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Abstract. In this review, I survey our current understanding of how the very 
first stars in the universe formed, with a focus on three main areas of interest: the 
formation of the first protogalaxies and the cooling of gas within them, the nature 
, and extent of fragmentation within the cool gas, and the physics - in particular the 

interplay between protostellar accretion and protostellar feedback - that serves to 
determine the final stellar mass. 
O |. In each of these areas, I have attempted to show how our thinking has developed 

D ■ over recent years, aided in large part by the increasing ease with which we can now 

' perform detailed numerical simulations of primordial star formation. I have also tried 

■ to indicate the areas where our understanding remains incomplete, and to identify 
' some of the most important unsolved problems. 
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^ . 1. Introduction 

For more than a decade - ever since the first release of the COBE results 
i-^ ■ (Mather et al., 1990; Smoot et al, 1992) - astrophysicists and cosmol- 

ogists have found themselves in the unusual situation of knowing more 
about the state of the universe when it was only 380,000 years old than 

■ when it was 200 million years old. Increasingly precise measurements 
of the cosmic microwave background (CMB), as exemplified by the 
recent results from WMAP (Bennett et al., 2003), together with a broad 
range of other observational constraints (Riess et al., 1998; Perlmutter 
et al., 1999; Bahcall et al, 1999; Valentine, Saunders, and Taylor, 



2000; Freedman et al., 2001; Percival et al, 2001; O'Meara et al., 
2001; Kirkman et al, 2003) have helped to confirm that we live in a flat 
universe, with approximately 5% of the closure density being provided 
by baryons, 25% by cold dark matter (CDM), and the remaining 70% 
by some form of 'dark energy' or cosmological constant. Models of such 
a universe - generally known as ACDM models - have been heavily 
studied for a number of years (see, e.g. Suginohara and Suto, 1991; 
Carroll, Press, and Turner, 1992; Gnedin, 1996ab) and many of their 
features are well understood. For instance, the evolution of the small 
inhomogeneities in the early universe that give rise to the observed 
temperature anisotropies in the CMB can be followed in great detail 
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(Seljak and Zaldarriaga, 1996), and the resulting predictions have been 
strongly confirmed by the WMAP results. 

The evolution of the dark matter component of the universe subse- 
quent to the epoch of last scattering at z ~ 1100 has also been studied 
intensively, using a wide range of techniques (see, for instance Seljak, 
2000; Benson et al, 2001; Cooray and Sheth, 2002). The general agree- 
ment between the results of these studies and an increasing number 
of observational tests (e.g. Gray et al, 2002) has lent further support 
to this overall picture, although some puzzles remain (Moore et al, 
1999; Navarro and Steinmetz, 2000). 

When it comes to understanding the behaviour of the baryonic 
component, however, we are on much shakier ground. Although cos- 
mological perturbation theory has given us a fairly good understanding 
of the behaviour of the baryons in the linear regime (Gnedin and Hui, 
1998; Meiksin, White, and Peacock, 1999; Singh and Ma, 2002), many 
details of the non-linear evolution of the baryons and the development 
of stars and galaxies are not understood. At the same time, we have 
little or no observational data to guide us. Although we now have obser- 
vational probes of the Universe at redshifts z > 6, thanks to the success 
of the Sloan Digital Sky Survey at finding high-redshift quasars (Fan 
et al, 2003), the strong metal lines observed in many of these quasars 
(Fan et al, 2001) are evidence that we are not yet probing the earliest 
epochs of star formation. 

Since observational limitations prevent us, for the time being, from 
directly studying the formation of the first stars and galaxies, work in 
this area has been primarily theoretical in nature. Although developing 
a theoretical understanding of primordial star formation may seem at 
first to be a hopelessly optimistic ambition - after all, there is still 
much that we do not understand about local star formation, despite the 
large quantity of observational data available to us - there are actually 
several good reasons to think that the problem may be a simpler one 
than understanding present day star formation. 

First, the initial conditions - small perturbations to a uniform cos- 
mological background - are simple and well understood (provided that 
the ACDM model remains as accurate on small scales as it has proved to 
be on the larger scales probed by galaxy surveys and the CMB). Second, 
the chemistry of primordial gas is also simple, at least in comparison 
to that of present day molecular clouds. It is therefore much easier to 
identify the critical reactions and to numerically simulate the chemical 
evolution of the gas. Third, magnetic fields, if present, are unlikely to be 
dynamically significant (Widrow, 2002); consequently, they are usually 
ignored. Finally, by restricting our attention to the first generation 
of stars to form, we can avoid the many complications posed by the 
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feedback of stars on their surroundings (see, for instance, Glover, 2001 
and references therein). 

Nevertheless, the problem remains a challenging one that involves 
processes occurring over a very wide range of length scales, from the cos- 
mological to the protostellar. A popular approach is to break this prob- 
lem up into a series of simpler problems, with different characteristic 
scales, that can be tackled individually. For instance: 

(i) When do the first protogalaxies 1 form, and how massive are they? 

(ii) How does gas evolve within these protogalaxies? Does it fragment, 
and if so, how large are the resulting fragments? When and why 
does fragmentation stop? 

(iii) What is the initial mass function (IMF) of the first stars, and what 
processes determine this? 

In this article, I review our progress at answering these questions. The 
body of the review is divided into three sections, each with a theme 
that broadly reflects one of the questions posed above, although there 
is inevitably a certain amount of overlap. 

There are, of course, many interesting questions concerning the first 
stars and galaxies which I have neither the time nor the space to 
properly address in this review; for instance, the question of how best 
to go about observing them; or the question of how they affect their 
environment both on small and large scales. Some of these questions 
are addressed in other recent reviews of early star formation (Barkana 
and Loeb, 2001; Bromm and Larson, 2004; Ciardi and Ferrara, 2004), 
which complement the material presented here. 

Throughout this paper, unless otherwise indicated, I adopt cosmo- 
logical parameters taken from the wmap concordance model (Spergel et 
ai, 2003). Specifically: Q m = 0.29, tt A = 0.71, Q h = 0.047, h = 0.72, 
ct 8 = 0.9, n s = 0.99. 

2. The formation of protogalaxies 

2.1. The first bound objects 

In CDM models, gravitationally bound objects form in a hierarchical, or 
'bottom-up' fashion, with the smallest, least massive objects forming 

1 A note on terminology: in this review, I use the term 'protogalaxy' as a con- 
venient shorthand for 'gravitationally bound gas cloud': the fact that something is 
described as a protogalaxy does not imply that it is actively forming stars, merely 
that it has the potential to do so. 
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first, and larger objects forming later through a mixture of mergers 
and accretion. The mass scale on which gravitationally bound objects 
begin to form (i.e. the minimum mass of a bound object) is set by the 
free streaming of the dark matter particles (Blumenthal et al., 1984). 
In general, this mass scale is many orders of magnitude smaller than 
scales of cosmological interest; for instance, in the neutralino model 
for CDM, M min ~ 1(T 7 M (Hofmann, Schwarz, and Stocker, 2001). 
However, the subsequent formation of larger objects occurs rapidly, 
and at most redshifts a large number of gravitationally bound objects 
(frequently referred to as 'dark matter halos') exist, with a wide range 
of different masses. 

Considerable effort has been devoted to determining the mass func- 
tion of dark matter halos as a function of redshift. The most widely 
used expression for the mass function is one originally suggested by 
Press and Schechter (1974): 

n(M, z) dM = \[^^4^7 exp ( ~—) dM. (1) 

Here n(M, z) dM is the comoving number density of halos at redshift z 
with dark matter masses in the interval (M, M + dM), p^m is the cos- 
mological background density of dark matter, and u = S c /[D(z)a(M)], 
where S c is a critical overdensity (generally taken to be 1.69), D{z) is the 
linear growth factor (Peebles, 1980; Carroll, Press, and Turner, 1992) 
and cr(M) is the rms fluctuation in the cosmological density field of 
dark matter smoothed on a mass scale M. A comprehensive discussion 
of the derivation of this equation is given in Bond et al. (1991). 

Comparisons with the results of N-body simulations at low redshift 
(Jenkins et al., 2001) and at high redshift (Jang-Condell and Hernquist, 
2001) demonstrate that the Press-Schechter mass function provides a 
reasonable fit to the true mass function, although a better fit to the 
simulation results can be obtained by using a modified form suggested 
by Sheth and Tormen (1999): 

„(M, ,) dM = A (l + -Ij) f-^ e,p (4) dM, (2) 

where v' = y/Ev, a = 0.707, A ~ 0.322 and q = 0.3. 

In either case, the basic form of the mass function is the same: it 
behaves as a power law for v <C 1 and falls off exponentially for v 3> 1. 
In CDM models, a(M) decreases monotonically with increasing mass, 
and so the most massive objects will also be the rarest. The transition 
to exponential behaviour occurs for v ~ 1, or a(M) ~ 5 C /D(z), and 
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so this transition occurs at a progressively smaller mass as we move to 
higher redshifts. 

Given a mass function of this type, is there any way to specify when 
the first halo of a given mass forms? Strictly speaking, the answer is 
no; the probability of finding a halo of any finite mass is never zero. 
In practice, however, we are more interested in determining when this 
probability grows to some interesting size, or when the number density 
of halos exceeds some specified threshold (which amounts to the same 
thing). This is most commonly done by specifying a value of v which is 
of interest; for instance, reference is often made to 3a halos, which are 
simply halos for which v = 3 and which therefore have a dark matter 
mass M satisfying: 

Such halos are moderately rare objects, representing no more than a few 
thousandths of the total cosmic mass (Mo and White, 2002), but are 
sufficiently common that one would expect to find many of them within 
a single Hubble volume. They are often taken to be representative of 
the earliest objects to form, although this choice is somewhat arbitrary. 

Unfortunately, while the Press-Schechter approach allows us to de- 
termine when the first dark matter halos of a given mass form, it 
does not, by itself, tell us when the first protogalaxies form, as it 
contains no information about the behaviour of the baryonic com- 
ponent of the universe. Unlike the dark matter, the baryons do not 
initially form structures on very small scales, since pressure forces act 
to suppress the growth of small-scale perturbations (Jeans, 1902; Jeans, 
1928; Bonnor, 1957). We can estimate the scale on which pressure forces 
become significant by equating the sound-crossing timescale, i sc , with 
the gravitational free-fall timescale, tg: if t sc < then perturbations 
can respond subsonically to changes in the gravitational field and will 
therefore remain in approximate hydrostatic equilibrium; on the other 
hand, if t sc > % then perturbations cannot respond subsonically, and 
some degree of gravitational collapse becomes inevitable. In gas with a 
density p and sound speed c s , we would therefore expect collapse to be 
suppressed on scales 

A more careful analysis using linear perturbation theory (Peebles, 1980) 
shows that in a purely baryonic universe, the growth of perturbations 
is completely suppressed on scales smaller than 
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where pb is the cosmological baryon density. This critical wavelength 
is commonly known as the Jeans length. The associated mass scale, 
known as the Jeans mass, is conventionally defined as 



3 



Mj = 3^HtJ • (6) 

The value of the Jeans mass depends on the baryon density, which is 
a simple function of redshift, and on the temperature of the intergalac- 
tic medium (through the dependence of Aj on c s ). The latter is simple 
to calculate at epochs prior to the onset of widespread star formation 
and is well approximated by (Galli and Palla, 1998) 

t = 410 (w) 2k < 7 > 

for redshifts z < 150. The corresponding Jeans mass at these redshifts 
is given by 

4.9 x 10 4 /l + z\ 3 / 2 w 

M ' = Ja^W bar) <8> 

To generalize this to the case of a universe containing both baryons 
and cold dark matter, it is tempting to simply replace the baryon 
density in the above equations with the total density p m = pb + Pdm, 
which would give us 

4.9 x 10 4 (\ + z\ 3 / 2 w 
M J= ,o rovTTo M « (9) 



(n m h 2 ) 1/2 V 150 

for z < 150; or in other words, a Jeans mass that is a factor (Qb/^m) 1 ^ 2 
smaller. In fact, the situation is not so simple, as perturbations can con- 
tinue to grow on small scales in the dark matter even when suppressed 
in the baryons. A linear treatment of this case is given in Gnedin and 
Hui (1998), but ultimately this treatment breaks down as small-scale 
structure in the dark matter begins to grow non-linearly. Although 
these non-linear effects have received little direct study, there is some 
evidence from numerical simulations that they can cause baryons to 
collapse on scales smaller than Aj (see, for instance, the discussion in 
section 2.1 of Haiman and Loeb, 1997), although any such collapse will 
be significantly delayed relative to the dark matter due to the influence 
of the gas pressure. In view of this, it is probably best to treat the 
value of Mj given by Equation (9) as an estimate of the scale on which 
pressure effects begin to dominate, rather than as an absolute lower 
limit to the protogalactic mass. 

Given Mj , we can go on to estimate the mass and formation redshift 
of the first protogalaxies by asking when the total mass of a 3a halo 
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Redshift (1+z) 

Figure 1. The evolution with redshift of M?, a (solid line), M^a (dot-dashed line) and 
Mj (dashed line). Protogalaxies will develop within 3<r dark matter halos once the 
mass of dark matter in the halo, Mz a , exceeds Mj; this occurs at z < 30. Similarly, 
protogalaxies will form in 4a halos once their dark matter mass, M^, exceeds Mj, 
which occurs at z < 40. 



first exceeds the Jeans mass. To properly answer this question, we would 
need to know the baryon fraction of these protogalaxies (i.e. their ratio 
of baryonic to dark matter). In practice, however, we know that this will 
be small and that the protogalactic mass will be dominated by the dark 
matter component. Therefore, for the purposes of a simple estimate it 
is sufficient to compare the Jeans mass with the dark matter mass of 
the 3a halo (hereafter M^ a ), which we can calculate using Equation (3). 
The evolution with redshift of both mass scales is plotted in figure 1. 
We can see from the figure that the first protogalaxies will have a total 
mass M ~ 10 4 M Q and will form at a redshift z ~ 30. It is also clear 
that uncertainties in Mj will have little effect on the estimated redshift, 
due to the sharp rise in M^ a with declining z. On the other hand, the 
use of a different criterion to identify our 'first' objects (e.g. considering 
4(j halos instead of 3a ones) has a rather larger effect on z, but has 
very little effect on the estimated protogalactic mass. 
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2.2. The importance of cooling 

Once a protogalaxy has formed, the next task is to determine how the 
gas within it evolves. In particular, we would like to know whether 
every protogalaxy that forms is capable of forming stars, or whether 
there are other prerequisites. 

We can gain considerable insight into this question by considering 
the thermal evolution of a parcel of gas that is undergoing gravita- 
tional collapse. The gravitational potential energy of the gas is trans- 
formed first into kinetic energy and thence into thermal energy through 
adiabatic compression, as well as the action of shocks if the flow is 
supersonic. Unless the gas can dissipate this thermal energy through ra- 
diative cooling, it must inevitably heat up. Since both density and tem- 
perature are rising, the pressure will increase rapidly, and ultimately 
will become large enough to halt the collapse. 

We can make this argument more quantitative by considering the 
gravitational stability of small perturbations within the collapsing gas. 
As in the cosmological case, we can derive a minimum unstable mass 
scale, again termed the Jeans mass, which scales as 

c 3 T 3 / 2 

m j « jj2 « -jn- ( 10 ) 

If we define an effective adiabatic index 

-, dlnT 

7efT = 1 + "Tj i 11 

din p 

then Mj will evolve with density as 

Mj oc p|(Teff-|)_ (12) 

Therefore, if 7 c ff > 4/3, the Jeans mass will increase during the collapse 
and will eventually become comparable to the mass of the protogalaxy, 
at which point collapse must halt. Since 7 er r = 5/3 for an atomic 
gas evolving adiabatically, it is clear that in the absence of radiative 
cooling, the increasing thermal pressure will bring collapse to an end 
long before protostellar densities are reached. Therefore, star formation 
is only possible if the gas can cool. 

The timescale on which cooling occurs is also of great importance. 
It has long been argued (Gott and Thuan, 1976; Rees and Ostriker, 
1977; Silk, 1977a) that the behaviour of gas in a collapsing protogalaxy 
depends upon the relative sizes of its cooling timescale, 

1 nkT 

" 7 _1A(T)' (13) 
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(where n is the particle number density and A(T) is the cooling rate 
per unit volume), its dynamical (or free-fall) timescale, given by 



td - = V^v (14) 

and the cosmological timescale, or Hubble time, 

( » = wry (15) 

It is easy to show that idyn is always less than tn, so there are only 
three possible arrangements: 

(i) *cool > *H > tdyn, 
(h) *H > icool > idyn, 
(hi) tn > idyn > *cool- 

In case (i), cooling takes place on a cosmological timescale, and is 
so slow that the gas evolves much as if there were no cooling at all. It 
quickly becomes pressure supported and remains so almost indefinitely, 
unless disturbed by an external event, such as a merger with another 
protogalaxy. In case (ii), the gas also becomes pressure supported, but 
subsequently contracts quasi-statically on a cosmologically interesting 
timescale. Finally, gas described by case (hi) never becomes pressure 
supported, but instead simply collapses at or near the free-fall rate. 

In practice, the situation is often far more complex than this analysis 
suggests, since the appropriate description for gas in a given proto- 
galaxy may vary with its location within the protogalaxy, and may 
also change over time as the density, temperature and/or chemical 
makeup of the gas change. Nevertheless, this scheme is a useful first 
approximation, and serves to further highlight the central role played 
by radiative cooling. 



2.3. Cooling and chemistry within primordial gas 



A number of potential cooling mechanisms exist within primordial gas 
(Anninos et ai, 1997), but many, such as Lyman-a cooling, operate 
only for T > 10 4 K, while the first protogalaxies have characteristic 
temperatures T ~ 100-1000 K. At these low temperatures, the domi- 
nant coolant is molecular hydrogen, H2, the most abundant primordial 
molecule. Therefore, to determine the cooling rate accurately, we re- 
quire an accurate value for the H2 abundance, which means that in 
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addition to studying the thermal evolution of the gas we must also 
study its chemical evolution. 

The chemistry of primordial gas has been investigated by a number 
of authors (Dalgarno and Lepp, 1987; Black, 1991; Abel et al, 1997; 
Galli and Palla, 1998; Standi, Lepp, and Dalgarno, 1998, 2002) and 
proves to be surprisingly complex despite the limited number of ele- 
ments involved. This complexity is due to the wide variety of different 
molecules and molecular ions that can be formed. However, if we are 
only interested in those aspects of the chemistry that affect the cooling 
rate, then we can make substantial simplifications (Abel et al., 1997): 
the chemical model can be reduced to a few processes that determine 
the ionization balance of the gas (e.g. collisional ionization, radiative 
recombination) , together with those reactions involved in the formation 
and destruction of H2 . 

The formation of H2 in local molecular clouds occurs primarily on 
the surface of interstellar dust grains: hydrogen atoms are adsorbed 
onto the surface of the grains, react to form H2 and subsequently 
escape back into the interstellar medium (Gould and Salpeter, 1963). 
In primordial gas, however, there is no dust, and so no possibility of 
forming H2 by this process. Instead, H2 formation is dominated by 
various sets of gas phase reactions. 

The simplest gas-phase reaction - direct radiative association of two 
hydrogen atoms to form H2: 



is strongly forbidden unless one of the hydrogen atoms is in an excited 
electronic state, and therefore plays an important role only in rather 
unusual circumstances, such as in the intergalactic medium near the end 
of the epoch of recombination (Latter and Black, 1991; Rawlings, Drew, 
and Barlow, 1993). It does not significantly influence protogalactic H2 
formation. 

Three-body formation of H2, via the reactions 



can play a significant role (Palla, Salpeter, and Stahler, 1983), but only 
at high densities (nu J> 10 8 cm -3 ), since the rate coefficients of these 
reactions are small. At lower densities, gas-phase formation of H2 is 
dominated by two sets of reactions. The first involves the H ion as an 
intermediate state 



H + H H 2 + 7, 



(Rl) 



H + H + H H 2 + H, 
H + H + H 2 H 2 + H 2 , 



(R2) 
(R3) 



H + e" -► H~ + 7, 
H +H H 2 + e", 



(R4) 
(R5) 
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and was first discussed in the context of the local ISM by McDowell 
(1961), and in a cosmological context by Peebles and Dicke (1968). The 
second set of reactions involves the H^" ion as an intermediary, and was 
first discussed in a cosmological context by Saslaw and Zipoy (1967) 



These two sets of reactions (hereafter the H pathway and the H^ 
pathway respectively) share two important characteristics. Firstly, both 
are limited by their initial step, since the radiative association reactions 
R4 and R6 occur at a much slower rate than the subsequent ion-neutral 
reactions R5 and R7. Secondly, the role played by free electrons in the 
H~ pathway is extremely similar to the role played by H + ions in 
the H2" pathway, and in both cases the H2 formation rate is directly 
proportional to the fractional ionization of the gas, provided that the 
latter is small. 2 

The main difference between the two pathways stems from the differ- 
ence in the rates of reactions R4 and R6: H~ forms via R4 much faster 
than forms via R6, and so the H~ pathway generally dominates the 
gas phase production of H2. 

The dependence of H2 formation on the presence of free electrons 
and protons might lead one to suppose that H2 formation will be very 
inefficient in low temperature gas, since the equilibrium ionization frac- 
tion is very low. In practice, however, moderate amounts of H2 can be 
formed if the protogalactic gas is not initially in ionization equilibrium. 
This is certainly the case in newly-formed protogalaxies, since the IGM 
itself is not in ionization equilibrium; instead, it retains a residual frac- 
tional ionization dating from the epoch of recombination. This comes 
about because the Hubble expansion ensures that the recombination 
timescale exceeds the expansion timescale before the IGM can reach 
equilibrium, freezing the fractional ionization at a value of approxi- 
mately 2 x 10 -4 (Stancil, Lepp, and Dalgarno, 1998). Protogalaxies 
forming from the IGM therefore begin with this small non-equilibrium 
fractional ionization. 

Once a protogalaxy has formed, this residual ionization quickly van- 
ishes, as the increased density leads to a greatly increased recombina- 
tion rate. However, there remains a brief window of opportunity in 
which H2 can form. Simple estimates of the resulting molecular frac- 
tion have been given by a number of authors (Susa et al, 1998; Nishi 

If the fractional ionization is large, then the mutual neutralization of H~ with 
H + and the dissociative recombination of become significant, and this simple 
relationship breaks down; this is discussed in more detail in Glover (2003). 



H + H+ - H+ + 7, 
H++H H 2 + H+. 



(R6) 
(R7) 
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and Susa, 1999; Oh and Haiman, 2002) and are typically in the range 
/h 2 = 10~ 3 -10~ 4 . For comparison, note that the molecular fraction in 
the IGM at this time is approximately 2 x 10~ 6 (Galli and Palla, 1998). 

Given the H2 abundance, density and temperature, it is then a 
simple matter to calculate the H2 cooling rate. Various parameteri- 
zations of this rate have been given in the literature; figure 2 shows 
some commonly cited examples, plotted as Ah 2 (^h™h 2 ) -1 ) where Ah 2 
is the H2 cooling rate per unit volume. The basic features of the cooling 
rate are straightforward: it falls off exponentially at low temperatures, 
due to the rather large excitation energy of the first accessible excited 
state (the J = 2 rotational state, which lies 512K above the J = 
para-hydrogen ground state), and is essentially negligible below 100 K; 
it scales with density as Ajj 2 oc n^ 2 at low densities, where radiative 
de-excitation dominates, and as Ah 2 oc «h 2 at high densities, where 
collisional de-excitation dominates and the level populations approach 
their local thermodynamic equilibrium (LTE) values. The transition 
between low density and high density behaviour occurs near a critical 
density n cr ~ 10 cm' 3 . 

The major uncertainty in the determination of the H2 cooling rate 
comes from uncertainties in the values of the collisional de-excitation 
rates, which are highly sensitive to the details of the potential energy 
surface used to calculate them, as well as to the method of calculation 
adopted (Lepp, Buch, and Dalgarno, 1995). As a result, there exists 
substantial disagreement in the literature on the form and magnitude 
of the H2 cooling rate at low densities, as can be seen from figure 2. 
However, recent calculations have removed much of this uncertainty, 
with the calculated collisional rates having more or less converged. 

2.4. Thermal evolution: simple models 

Armed with an appropriate set of chemical reaction rates and an ac- 
curate H2 cooling rate (see, for example, Abel et al, 1997 or Glover, 
2001), the next step is to follow the coupled chemical, thermal and 
dynamical evolution of a protogalaxy as it forms in order to determine 
its fate. 

The simplest approach to this problem dispenses entirely with any 
attempt to accurately simulate the dynamical evolution of the proto- 
galaxy. Instead, the density evolution is specified in advance, and the 
model focuses on determining the chemical and thermal evolution of 
the gas. For instance, it is frequently assumed that if cooling is effective, 
then the density evolution will be the same as in pressure-free collapse. 
This approximation relies on the assumption that pressure gradients 
are everywhere small compared to gravitational forces. 
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Temperature (K) 

Figure 2. A comparison of various parameterizations of the H2 cooling function, 
plotted in units of erg cm 3 s _1 . Rates are computed assuming that nu S> iih 2 , and 
that the ortho-to-para ratio is 3:1. The lower set of lines corresponds to a gas density 
n,H = 10 6 cm -3 ; the upper set corresponds to iin — 10° cm~ 3 . Solid line - Le Bourlot, 
Pineau des Forets, and Flower (1999); dashed line - Galli and Palla (1998); dotted 
line - Lepp and Shull (1983); dash-dotted line - Hollenbach and McKee (1979). 



A number of authors have considered the problem of protogalac- 
tic collapse within this framework (Matsuda, Sato, and Takeda, 1969; 
Hutchins, 1976; Yoshii and Sabano 1979, 1980; Carlberg, 1981; Palla, 
Salpeter, and Stahler, 1983; Villere and Bodenheimer, 1987; Susa, Ue- 
hara, and Nishi, 1996; Omukai, 2000; Flower and Pineau des Forets, 
2001), often supplementing it with the additional assumptions of spheri- 
cal symmetry and uniform density Use of these approximations reduces 
the problem to one of computing the chemical and thermal evolution 
of a single representative parcel of gas. 

Most of these models predict the same general type of behaviour. 
Initially, the H2 cooling rate is negligible and the evolution of the gas is 
very close to adiabatic. As the collapse proceeds, however, the increas- 
ing temperature, density and H2 abundance all combine to dramatically 
increase the H2 cooling rate and decrease t coo i- Eventually, t coo \ becomes 
comparable to the collapse timescale, and the collapse ceases to be 
even approximately adiabatic. Instead, the temperature reaches a peak 
and then decreases at higher densities as radiative cooling becomes 
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increasingly dominant over compressional heating. The quantitative 
details, such as the value of the peak temperature, are sensitive to 
the treatment of H2 cooling and gas chemistry adopted, and generally 
vary from model to model, although never by more than a factor of a 
few. 

The only case in which this type of model predicts substantially 
different behaviour is when some other process, such as UV photodis- 
sociation, acts to reduce the H2 abundance (see, e.g. Omukai, 2001). 
In this case, the protogalaxy may be unable to form sufficient H2 to 
cool the gas before it reaches a temperature and density at which 
collisional dissociation of H2 becomes significant. This results in the 
gas temperature continuing to rise until the onset of Lyman-a cooling 
at a temperature of approximately 10 4 K. 

An alternative model is presented by Tegmark et al. (1997). They 
make a similar set of approximations (spherical symmetry, uniform den- 
sity, free-fall collapse) , but halt the collapse when one of two conditions 
is met: 

(i) The gas temperature exceeds the virial temperature of the proto- 
galaxy, defined as (Blanchard, Valls-Gabaud, and Mamon, 1992): 

Tvir = GMiimy 

2 k -f^vir 

(ii) The mean density of the gas exceeds the mean density of the dark 
matter halo. The latter can be written as 



where A = 187r 2 for an Einstein-de Sitter cosmology (Peebles, 
1980); analogous values for open or A-dominated cosmological mod- 
els are given in Bryan and Norman (1998). 

If condition (ii) is met, then shocks are assumed to raise the gas 
temperature instantaneously to T vir at the end of the collapse. 

Tegmark et al. (1997) make no attempt to follow the further dynam- 
ical evolution of the gas. Instead, they hold its density constant, and 
study its subsequent thermal and chemical evolution. If the gas tem- 
perature decreases by more than 25% during an interval corresponding 
to a 25% decrease in redshift, i.e. if 



where z c is the redshift at which the collapse terminates, then the pro- 
togalaxy is considered to be able to cool effectively. Otherwise, cooling 
is considered to be ineffective. 



PDM = (1 + A)p dm , 



(17) 



T(0.75z c ) < 0.75T(z c ), 



(18) 
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By performing this analysis for protogalaxies with a wide range of 
masses and collapse redshifts, Tegmark et al. are able to map out the 
region in M-z c parameter space corresponding to protogalaxies that 
can cool effectively. They find that at each redshift they can identify 
a minimum mass M imn such that protogalaxies with M > M m j n cool 
effectively, while those with M < M m \ n do not. For instance, at z = 30, 
they find that M nun ~ 10 6 M Q , two orders of magnitude larger than 
the Jeans mass at that redshift. This result can also be expressed in 
terms of a minimum virial temperature, T m i n , related to M m i n through 
Equation (16); at z = 30, this is 1000 K. While the values of M m i n and 
Tmin obtained in this way are somewhat sensitive to the choice of H2 
cooling function (Abel et al, 1998; Glover, 2001), the basic behaviour 
remains the same. 

These two approaches - the free-fall collapse model and the Tegmark 
et al. model - therefore present us with two distinct scenarios for the 
formation of protogalaxies. The free-fall models predict that every pro- 
togalaxy can cool, with the onset of cooling occurring once the gas 
has reached a temperature of approximately 1000 K (give or take a 
factor of two). On the other hand, the Tegmark et al. model predicts 
that only those protogalaxies with T v ; r > 1000 K will cool; smaller 
protogalaxies, with lower virial temperatures, will simply remain as 
pressure-supported gas clouds and will not form stars. 

Are either of these scenarios correct? Ultimately, this depends upon 
whether the approximations on which they are based are justified. This 
is a question that is best addressed through the use of more detailed 
numerical simulations. 

2.5. Thermal evolution: numerical simulations 

The biggest problem that we face when trying to simulate protogalactic 
collapse numerically is the wide range of length scales that we are 
required to resolve. For example, consider the collapse of a 10 6 M 
protogalaxy. This has a characteristic size (as given by the virial radius) 
of approximately 3 kpc in comoving units, corresponding to 100 pc in 
physical units at z = 30. To properly simulate its cosmological envi- 
ronment, we should follow the evolution of the gas and dark matter 
on scales that are two to three orders of magnitude larger (see, for 
instance, the resolution study of Ricotti, Gnedin, and Shull, 2002a), 
while to resolve star formation within it, we need to be able to follow 
the gas down to scales of the order of an AU. The total dynamical 
range required in order to resolve all of this within a single simulation 
is therefore approximately 10 10 , many orders of magnitude larger than 
can be covered with a single fixed grid. 
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The simplest way in which we can obtain the required dynamical 
range is to use a Lagrangian grid, i.e. one which moves with the fluid 
flow. This is particularly effective if one assumes that protogalaxies 
are spherically symmetric, as in this case we can use a simple one- 
dimensional Lagrangian code, such as that described by Thoul and 
Weinberg (1995). 

The earliest studies of this type were performed by Peebles and Dicke 
(1968) and Poppel (1975), but the initial conditions that they adopted 
- isolated, isothermal clouds, initially in hydrostatic equilibrium - are 
not appropriate for protogalaxies forming by dynamical collapse within 
an expanding cosmological model. More recent simulations by Haiman, 
Thoul, and Loeb (1996), Oliveira et al. (1998ab) and Stachniewicz 
and Kutschera (2003) begin from more appropriate cosmological initial 
conditions and all three groups obtain broadly similar results. 

The most significant result of these simulations is the demonstration 
that the dynamical evolution of a small, H2 cooled protogalaxy is not 
well described by pressure-free gravitational collapse. Instead, pressure 
plays a important role, particularly for protogalaxies with masses near 
Mj. It has two main effects. In the initial stages of collapse, when the 
flow is subsonic, it delays the collapse of the gas relative to the dark 
matter, resulting in a density profile which is less centrally concentrated 
than would otherwise be the case. At a later time, once the infall has 
become supersonic, the finite pressure leads to the formation of an 
accretion shock near the virial radius of the halo. The majority of the 
H2 that forms does so in the post-shock gas, and it is the conditions 
there that determine whether or not the protogalaxy is able to cool 
effectively (Haiman, Thoul, and Loeb, 1996). 

Of the two scenarios discussed in the previous section, the Tegmark 
et al. model clearly provides the better description. The major point of 
disagreement is the protogalactic density profile: for simplicity, Tegmark 
et al. assume a uniform density profile, while the simulation results 
show that the true profile is much closer to that of an isothermal sphere. 

Unfortunately, one-dimensional Lagrangian simulations, although 
simple to perform and quick to run, ultimately give us a rather lim- 
ited view of protogalactic formation, since they do not include many 
important physical effects such as rotation and turbulence. To identify 
the role that these factors play, we need to use fully three-dimensional 
hydrodynamical simulations. 

The move to three dimensions necessitates a change in our compu- 
tational strategy, as grid-based Lagrangian codes do not handle three 
dimensional flows well due to the severe grid distortion that tends to 
occur and which causes a dramatic loss of accuracy. This problem can 
be circumvented to some degree through the use of 'hybrid' codes, 
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which switch to Eulerian (i.e. fixed) coordinates in regions of high 
deformation. Examples include the codes of Gnedin (1995; see also 
Gnedin and Bertschinger, 1996) and Pen (1998). However, although 
this technique has been used to study protogalactic feedback (Ostriker 
and Gnedin, 1996; Gnedin and Ostriker, 1997; Ricotti, Gnedin, and 
Shull, 2002ab), it has not been used to study primordial star formation 
in any detail. 

Another way to avoid the grid distortion problem is to abandon 
the use of a grid, and to switch to a particle-based Lagrangian tech- 
nique such as smoothed particle hydrodynamics (SPH; see Gingold 
and Monaghan, 1977; Lucy, 1977; Monaghan, 1992). Alternatively, the 
required dynamical range can be obtained using fixed grids if multiple 
nested grids, or some form of grid refinement are used. Both of these 
approaches are discussed in more detail below. 



2.5.1. SPH simulations 

Several authors have studied the formation of protogalaxies using SPH 
(Bromm, Coppi, and Larson, 1999, 2002; Fuller and Couchman, 2000; 
Yoshida et al, 2003). Fuller and Couchman (2000) used the hydra 
cosmological SPH code (Couchman, Thomas, and Pearce, 1995) to 
study the collapse of uniform, spherical protogalaxies of various masses 
at a range of different redshifts, in order to test the predictions of 
Tegmark et al. (1997). They obtained broadly similar results, although 
their values of T m [ n are roughly a factor of two smaller than those of 
Tegmark et al., a discrepancy which may simply be due to the different 
H2 cooling functions used in the two studies. 

Fuller and Couchman also studied protogalactic formation in a more 
realistic cosmological simulation. They demonstrated that the evolution 
of the most massive object in their simulations could be divided into 
two main phases. In the first phase, the halo mass is less than the Jeans 
mass, pressure forces dominate, and the baryonic overdensity is small, 
but non-zero. This phase corresponds to the delayed collapse phase seen 
in the one-dimensional simulations discussed above. 

During this phase, the mass of the halo continues to increase, driven 
to a large extent by mergers with smaller dark matter halos. As the 
mass nears Mj, the gas density profile begins to steepen significantly as 
gravitational forces become dominant. At this stage, the gas tempera- 
ture is already significantly higher than the temperature of the IGM, 
while the fractional H2 abundance in the central dense region is of order 
10~ 4 , two orders of magnitude larger than its initial value. Nevertheless, 
the cooling time remains longer than the dynamical timescale, and the 
gas is not yet self-gravitating. 
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The second phase begins once t coo i drops below idyn at the centre of 
the protogalaxy. The central gas rapidly cools to T ~ 150 K, and the 
consequent reduction in pressure support leads to a substantial increase 
in the central density. The gas eventually becomes self-gravitating at 
z ~ 20, at which point the protogalaxy has a mass M ~ 4 x 10 5 M Q , 
close to the estimate of M m \ a at that redshift. 

Inspection of the spherically averaged temperature and density pro- 
files of the protogalaxy allows us to identify several distinct regions. 
From the outside in, we have: 

(i) Cosmological infall, terminated by an accretion shock 

(ii) A broad, post-shock region, where heating from adiabatic compres- 
sion competes with H2 cooling, and where t coo i > tdyn 

(iii) A cold, dense central core, where i coo i < tdyn- 

Unfortunately, the structure of the core region is not well resolved in 
Fuller and Couchman's simulation, since its extent is comparable to 
the minimum smoothing length of their SPH code. 

Bromm, Coppi, and Larson (1999, 2002) simulated protogalactic 
collapse using a modified version of the treesph code (Hernquist and 
Katz, 1989). They concentrated on following in detail the evolution of 
a single protogalaxy and consequently adopted simplified initial con- 
ditions: a spherical overdensity, set into rigid rotation with a specified 
angular momentum and perturbed on small scales using the Zeldovich 
approximation (Zeldovich, 1970) with a P{k) oc fc~ 3 power spectrum. 
By focusing on a single protogalaxy and neglecting its cosmological 
environment, they were able to follow its collapse to high densities 
[n < 10 8 cm" 3 ). On large scales, their results confirm those of Fuller 
and Couchman (2000): the gas initially evolves adiabatically, is heated 
up to T ^ T v [^ in an accretion shock, and subsequently cools to T ^ 
100-200 K in the dense central regions. On small scales, their greater 
resolution allowed them to follow the formation of structure within the 
dense, cold gas. This portion of their results is discussed in detail in 
section 3. 

Finally, Yoshida et al. (2003) used the gadget code (Springel, 
Yoshida, and White, 2001) to explore the cosmological environment 
in which the first protogalaxies form. They performed the largest pro- 
togalactic simulation to date, using 48 million SPH particles to simu- 
late the evolution of a cosmological volume which is 6OO/1" 1 comoving 
kiloparsecs on a side. This simulation had a mass resolution of approxi- 
mately 5000ft," 1 M and a spatial resolution of approximately 50/i 1 pc, 
and so resolved little of the internal structure of the protogalaxies. 
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On the other hand, it did allow a large sample of protogalaxies to be 
studied within a single consistent simulation, and was therefore a useful 
complement to more detailed studies of single protogalaxies. 

Yoshida et al. found that to cool efficiently, protogalaxies in their 
simulation must have masses M ~ 5 x lO 5 /},^ 1 M Q , and fractional H2 
abundances /h 2 — 2 x 10~ 4 . Moreover, while all of the protogalaxies 
with large H2 abundances also had large masses, the converse was not 
true; some protogalaxies with masses above 5 x 10 5 /i _1 M Q did not 
form enough H2 to cool. Further investigation of these protogalaxies 
showed that they were gaining mass more rapidly than their cooler 
counterparts, leading Yoshida et al. to suggest that their temperatures 
were being kept high by the compressional heating associated with 
frequent merger activity. If this is the case, then it implies that the 
ability of a particular protogalaxy to cool and form stars depends on 
its dynamical history as well as its current mass. Further investigation 
of this point is clearly warranted. 



2.5.2. Multigrid simulations 

The basic idea behind a multiple grid (or multigrid) Eulerian simulation 
is to take a single top-level grid that is large enough to represent the 
whole volume of interest, and then to supplement it with one or more 
levels of subgrids in regions in which higher resolution is desired. Since 
much of the volume of a typical cosmological simulation is filled with 
under-dense material that is well resolved by the top-level grid alone, 
this technique can dramatically improve the dynamical range of an 
Eulerian simulation for only a small increase in its computational cost. 

In the simplest implementation of the multigrid technique, the place- 
ment of the grids is specified at the beginning of the simulation and does 
not subsequently alter. This is the approach used in the protogalactic 
simulations of Abel et al. (1998). They used the hercules code (An- 
ninos, Norman, and Clarke, 1994; Anninos et al., 1997) to study the 
growth of 3a and 4a density peaks within a large cosmological volume. 
The simulations were performed using a top-level grid with a resolution 
of 128 3 , together with a quarter-size subgrid with the same resolution, 
for an effective resolution of 512 3 . The initial conditions were arranged 
to ensure that the density peak would remain within the region covered 
by the subgrid. 

These simulations were able to resolve the basic filamentary struc- 
ture of the IGM surrounding the protogalaxies, and gave some indica- 
tions that gas cooling within the protogalaxies was not particularly ef- 
ficient. However, the protogalaxies remained under-resolved, rendering 
these conclusions uncertain. 
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A natural way to improve the resolution would be to add more 
subgrids. However, as one does this, it becomes increasingly difficult to 
ensure that the subgrids are placed correctly, since at the beginning of 
the simulation it is generally not possible to determine exactly which 
regions will require very high resolution. Fortunately, this problem can 
be overcome by the use of a technique called adaptive mesh refinement. 

In an adaptive mesh simulation, the placement of subgrids is not 
specified a priori. Instead, one or more refinement criteria are specified, 
and local subgrids are created as required to ensure that these criteria 
are always satisfied. Adaptive mesh codes have been used with great 
success in a number of areas of astrophysics, as discussed in the recent 
review of Norman (2004). Their use in the study of primordial star 
formation was pioneered by Abel, Bryan, and Norman (2000, 2002) in 
a pair of highly influential papers. 

In the first of these papers, Abel, Bryan, and Norman (2000) used 
the ENZO code of Bryan and Norman (1997ab) to follow the evolution 
of protogalactic gas from cosmological scales down to densities of order 
10 6 cm -3 . They achieved a maximum resolution of 0.4 pc (in comoving 
units) within a box that was 128 comoving kpc on a side, for a total 
dynamic range of approximately 3 x 10 5 . This simulation produced a 
protogalaxy with the same basic temperature profile as found in the 
SPH simulations: a cold infall region, an accretion shock at r ~ r vrr , 
a subsequent broad cooling zone, and a cold, dense central region. 
In a subsequent simulation, described in Abel, Bryan, and Norman 
(2002), they included additional molecular physics (the three-body for- 
mation of H2) and followed the collapse to significantly higher densities, 
eventually reaching a minimum physical resolution of a few tens of AU. 

Although Abel, Bryan, and Norman present many of their results, 
such as the temperature and density profiles, in the form of spherically 
averaged quantities, they also present a number of slices through their 
simulations on different scales. These demonstrate that the assumption 
of spherical symmetry is a relatively crude approximation, particularly 
on large scales, since much of the gas falling into the potential well of 
the protogalaxy does so along a few overdense filaments, rather than 
in a spherically symmetric fashion. 

The very high dynamical range achieved in their simulations also 
allowed Abel, Bryan, and Norman to study in detail the evolution of 
the dense gas at the centre of the simulated protogalaxy. This portion 
of their results is discussed later, in section 3. 
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2.6. Summary 

By combining the detailed results of the numerical simulations de- 
scribed in the previous section with the more general theoretical con- 
siderations of sections 2.1-2.4, we are able to put together a reasonably 
comprehensive answer to the first of the questions posed in the in- 
troduction: when do the first protogalaxies form, and how large are 
they? 

As we have seen, the earliest protogalaxies will form at a redshift 
of 30-40, and will have masses of order 10 4 M©. However, these proto- 
galaxies will form little H2 and will not be able to cool effectively. They 
are therefore extremely unlikely to form stars. The earliest star-forming 
protogalaxies will form later, at z ~ 30, and will be more massive, with 
masses of order 10 5 -10 6 M©, and virial temperatures of order 1000 K. 

Finally, it should be noted that these results assume a CDM-based 
cosmological model. If this turns out to be an incorrect description of 
dark matter on small scales, then we should expect these numbers to 
change significantly. For instance, if some form of warm dark matter is a 
more appropriate description, then the first protogalaxies will be larger 
(M ~ 10 7 M Q ) and will form at lower redshift (z ~ 20), as demonstrated 
in the recent simulation by Yoshida, Sokasian, Hernquist, and Springel 
(2003). However, models such as this have great difficulty accounting 
for the high electron scattering optical depth detected by WMAP (Kogut 
et al, 2003), and at present there seems little reason to prefer them 
over CDM. 



3. Fragmentation 

Once we have established that a significant amount of protogalactic 
gas can cool and condense on a cosmologically interesting timescale, 
the next step is to investigate what happens to this gas. In particular, 
we would like to know whether any of it forms stars, and, if so, how 
many stars form, over what timescale, and with what IMF? 

Crucial to determining this is an understanding of the degree to 
which the protogalactic gas fragments during its dynamical evolution: 
does the cold gas sink into the centre of the halo, forming a single mas- 
sive clump? Or does it break up into many smaller clumps? To address 
these questions, in the following sections I examine the effectiveness 
of the various mechanisms that may bring about fragmentation, and 
discuss the results of the most recent numerical investigations. 
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3.1. Hierarchical fragmentation 



An obvious place to begin is with the force responsible for the formation 
of the protogalaxy itself: gravity. For gravitational fragmentation to be 
effective, two important conditions must be met. First, any fragments 
that form must be gravitationally bound. Second, regions that begin 
to collapse due to their own self-gravity must be able to complete this 
collapse on a timescale shorter than the dynamical timescale of the flow 
in which they are embedded; otherwise, they will be disrupted before 
they have time to grow into distinct objects. 

The question of whether a particular fragment is gravitationally 
bound depends, in the general case, on a number of factors: the mass of 
the fragment, its internal velocity field and pressure distribution, the 
properties of the surrounding gas flow etc. (see McKee and Zweibel, 
1992 for a detailed discussion). However, much of the work done on 
gravitational fragmentation in a primordial context makes the simpli- 
fying assumptions that the only forces acting are gravity and pressure, 
and that the latter can be neglected on scales larger than the Jeans 
length. Although the validity of these assumptions is questionable, they 
provide a simple starting point for an investigation of protogalactic 
fragmentation, so I will briefly discuss the conclusions they lead us to, 
before going on to consider more detailed models. 

My starting point is the work of Hoyle (1953). He considered the 
gravitational collapse of a homogeneous protogalaxy, and showed that 
on scales where pressure can be neglected, the second of our condi- 
tions for fragmentation will always be satisfied. His argument is very 
simple: gas with a density p collapses gravitationally on a timescale 
iff oc (G/?) -1 / 2 , while a perturbed region with a density p' will collapse 
on a timescale t' s oc (Gp') -1 / 2 . If p' > p, then t' s < %, and so the 
perturbed region will collapse faster than the main flow. In the proto- 
galactic case, this implies that overdense regions within the protogalaxy 
will be able to collapse under their own self-gravity in less time than it 
takes for the protogalaxy itself to collapse. Therefore, the protogalaxy 
will fragment. 

Hoyle also pointed out that one can apply precisely the same argu- 
ment to every fragment that forms within the protogalaxy: provided 
that they contains overdense regions, and that the neglect of pres- 
sure forces remains appropriate, they too will fragment (as will the 
fragments of these fragments etc.). Hoyle therefore argued that the pro- 
togalactic gas would continue to fragment on smaller and smaller scales 
until pressure forces finally intervened to prevent further fragmentation, 
a scenario that has come to be known as hierarchical fragmentation. 
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Hoyle's original semi-quantitative argument was subsequently placed 
on a sounder mathematical basis by Hunter (1962), who performed 
a linear perturbation analysis of a uniform, spherical, pressure-free 
collapse and showed that any overdense perturbation would be gravita- 
tionally unstable and would quickly grow until it became nonlinear in 
less than a free-fall time. Hunter (1964) expanded on this analysis by 
considering second order terms and showed that these would accelerate 
the growth of overdense regions. Similar analyses have also been made 
by Savedoff and Vila (1962) and Silk (1982), with similar results. 



3.2. The opacity limit 



An important prediction of the hierarchical fragmentation scenario is 
that the smallest fragments will have sizes of the order of the Jeans 
length, and hence masses of the order of the Jeans mass, since this is 
the scale on which pressure balances gravity However, both Aj and 
Mj are functions of the density and temperature of the gas, and will 
change as the protogalaxy evolves. To identify the minimum fragment 
mass, we must therefore determine how small Mj becomes during the 
collapse. 

Recall that we can write Mj in terms of the gas density as 

Mj oc pi^ff-D, (19) 

where the effective adiabatic index ^y e g is 

dlnT . . 

7erT = 1 + ~Tj • (20) 

dm p 

For 7 c ff < 4/3, the Jeans mass decreases with increasing density, while 
for 7 c fj > 4/3 it increases. Hoyle (1953) suggested that a transition 
from 7 c ff < 4/3 to 7 c fr > 4/3 would occur when the gas first became 
optically thick, under the assumption that this would mark a change 
from isothermal evolution to adiabatic evolution. The density and tem- 
perature of the gas at this time would then set the minimum fragment 
mass. This basic idea - that it is the opacity of the gas which sets a 
lower limit on the mass of a fragment and therefore on the mass of 
a star - has become known as opacity-limited fragmentation, and has 
been studied by a number of authors. 

Low and Lynden-Bell (1976) and Silk (1977b) both follow Hoyle 
in assuming that the minimum mass is reached at the moment that 
a fragment first becomes optically thick. They also assume that the 
fragment is in thermal balance at this time, with compressional heating 
balanced by radiative cooling. These assumptions provide them with 
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two equations relating T F , p F and k f (the temperature, density and 
opacity of the fragment at the moment that it becomes optically thick): 

^j(Pf,T f ) 

k fPf 2 = ' ^ ' 

T C ( PF ,T F ) = A r {p F ,T F , KF ), (22) 

where T c is the compressional heating rate and A r is the radiative 
cooling rate. 

We can use these equations to express the minimum fragment mass 
M F in terms of a single unknown - for instance, Low and Lynden-Bell 
(1976) write it in terms of the opacity as 

M F = 2.5 x 10- V 16/T ( 1A M , (23) 

\K F J 

where p is the mean molecular weight and kq is the opacity due to 
Thomson scattering - but to fully determine M F , we need an additional 
piece of information. In Hoyle's original analysis, this comes from the 
assumption that the radiative cooling is dominated by Lyman-a emis- 
sion, which fixes the temperature at approximately 10 4 K. On the other 
hand, Silk (1977b) considers a case where the cooling and opacity are 
both dominated by dust, in which case the observed dust temperature 
provides the additional information required. In general, however, we 
can only determine M F if we know something of the previous thermal 
history of the gas. 

Rees (1976) studied the opacity limit from a different viewpoint, 
arguing that the essential requirement for continued isothermal evolu- 
tion is that a fragment be able to radiate away its gravitational binding 
energy in less than a free-fall time, and that opacity is only important 
inasmuch as it limits the maximum radiative rate, which cannot exceed 
that of a black body of a similar temperature. This fact can be used to 
derive the minimum temperature that a fragment must have in order 
to radiate sufficient energy, which in turn can be used to determine 
M F : 

M F = M c p-^r 1/2 (—X\ (24) 

where M c is the Chandrasekhar mass, m p is the mass of a proton 
and / is a radiative efficiency factor, defined as the ratio of the actual 
radiation rate to the black-body rate: 

_ jF v dv 
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From Equation (24), we see that if / ~ 1, then Mp will be of the order 
of a solar mass (or less, if Tp is very small), while if / <C 1, as may 
occur if the cooling is dominated by a few narrow emission lines, then 
Mp may be of the order of tens or hundreds of solar masses. 

More recently, Masunaga and Inutsuka (1999) have reconsidered the 
conditions under which isothermal evolution comes to an end. They 
show that this will inevitably occur once the radiative cooling rate 
becomes unable to keep pace with the compressional heating rate, and 
that this may take place in either the optically thin or optically thick 
regime, depending on the details of the collapse, but is unlikely to 
coincide with the instant at which r = 1. This suggests that a better 
procedure for determining Mp is to follow the actual thermal evolution 
of the gas. 



3.3. Simple numerical models 



Various authors have attempted to calculate Mp by modelling the ther- 
mal evolution of the collapsing protogalactic gas. One of the earliest 
was Yoneyama (1972), who constructed an evolutionary track for the 
gas in density-temperature space by assuming that it always satisfied 
the following conditions: 

icool = (26) 
t H2 = max(t ff ,t rcc ), (27) 

where tn 2 is the H2 formation timescale, given by 

tn 2 = r " Ha , (28) 
k H - n c n H 

where k H - is the rate coefficient for the formation of H~ by radiative 
association of electrons and Hi (reaction R4), and t rec is the recombi- 
nation timescale, given by 

^rec T j (29) 

n-rec 

where k rec is the rate coefficient for radiative recombination. Yoneyama 
used this technique to calculate the evolution of Mj until either the gas 
became optically thick or became hot enough to collisionally dissociate 
H2. Mp could then be computed simply by finding the minimum value 
of Mj reached along the evolutionary trajectory. Yoneyama found that 
in small protogalaxies, the minimum value was reached shortly before 
the gas became hot enough to collisionally dissociate H2 and that M-p ~ 
60 M ; in larger protogalaxies, the greater column density of H2 caused 
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the gas to become optically thick before collisional dissociation could 
occur, and the resulting value of Mp was considerably larger. 

A more common approach is to specify the form of p(t) in advance, 
often by constructing some extremely simplified dynamical model for 
the protogalaxy, and then to use this as an input into a more detailed 
chemical and thermal model. Many such models exist (Hutchins, 1976; 
Silk, 1977a; Carlberg, 1981; Hasegawa, Yoshii, and Sabano, 1981; Palla, 
Salpeter, and Stahler, 1983; Lepp and Shull, 1984; Lahav, 1986; Villere 
and Bodenheimer, 1987; de Araujo and Opher, 1989; Susa, Uehara, and 
Nishi, 1996); I will discuss only a few notable examples. 

Hutchins (1976) studied the thermal evolution of a variety of spheri- 
cal and spheroidal protogalaxies using a very simplified chemical model 
consisting of only four reactions - formation of H2 via the H~ pathway 
(reactions R4-R5), plus radiative recombination of hydrogen 

H+ + e" H + 7, (R8) 

and collisional dissociation of H2 by H: 

H 2 + H -> 3H. (R9) 

His cooling function was equally simple, and included only Compton 
cooling and H2 rotational cooling. The calculations were terminated 
once H2 dissociated. Hutchins found a minimum fragment mass of 
approximately 200 Mq, within a factor of a few of Yoneyama's result. 

Silk (1977a) and Carlberg (1981) both investigated the effects of 
including Lyman-a cooling in the model, and showed that it has a 
dramatic effect. The reason is that the collisional dissociation of H2 no 
longer results in a permanent transition to adiabatic evolution. Instead, 
the gas heats up adiabatically for a short while, until its temperature 
reaches 10 4 K, following which it again begins to evolve isothermally. 
This second phase of isothermal evolution is eventually terminated once 
the fragments become optically thick in the continuum. It allows Mp 
to reach much smaller values than would otherwise be possible. Silk 
(1977a) estimated the minimum fragment mass to be approximately 
0.3 M ; Carlberg (1981), using a more detailed treatment, found Mp ~ 
0.5 M . 

Finally, Palla, Salpeter, and Stahler (1983) investigated the impor- 
tance of three-body H2 formation (reactions R2-R3) and showed that 
at high densities (n > 10 8 cm -3 ), these reactions become very effec- 
tive, rapidly converting the bulk of the hydrogen to molecular form. 
This dramatically increases the H2 cooling rate, delaying the collisional 
dissociation of H2 until much higher densities are reached, and conse- 
quently lowers Mp. For the representative example of a 5 x 10 4 M 
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cloud, Palla, Salpeter, and Stahler (1983) found a minimum fragment 
mass of Mp ~ 0.1 M Q . 

The main lesson to draw from these attempts is the importance of an 
accurate treatment of the thermal evolution of the gas, which requires 
a comprehensive treatment of the microphysics - comparison of the 
results of Hutchins (1976) with those of Silk (1977a) or Palla, Salpeter, 
and Stahler (1983) demonstrates the danger of using an oversimplified 
chemical or thermal model. However, the specific predictions of these 
models rely upon the accuracy of the assumptions made regarding the 
density evolution, and we have already seen that simple collapse models 
generally do not perform well in this respect. More importantly, these 
predictions rely on the correctness of the basic assumptions underlying 
the hierarchical fragmentation scenario. Unfortunately, there are good 
reasons to believe that these assumptions are not correct, as I discuss 
in the next section. 



3.4. The case against hierarchical fragmentation 



Recall that the hierarchical fragmentation scenario is based on two 
major assumptions: first, that the balance between gravity and pres- 
sure is the only determinant of whether a fragment is gravitationally 
bound, and second, that the gas flow on scales larger than the Jeans 
length can be approximated by pure gravitational free-fall. The first 
of these assumptions implies that all fragments with M > Mj will be 
gravitationally bound; the second, that the gas can quickly fragment 
down to the Jeans scale, provided that it starts from moderately uni- 
form initial conditions. Clearly, both of these assumptions represent 
substantial simplifications. However, the real question is whether the 
simplified picture remains accurate; in other words, did hierarchical 
fragmentation actually occur in real protogalaxies? 

One important argument against this scenario was first advanced by 
Layzer (1963). He argued that as the protogalactic gas fragmented, the 
individual fragments would tend to acquire angular momentum from 
the gravitational torques exerted on them by their neighbours. If this 
angular momentum was subsequently conserved during the evolution of 
the fragment, then it would limit the extent to which it could contract, 
and would help to stabilize it against further fragmentation. 

A convenient way to parameterize this is in terms of the dimension- 
less spin parameter 

3\E\^ 2 v rot 

a =gm^"^' (30) 

where J and E represent the fragment's total angular momentum and 
total energy respectively. If angular momentum is conserved, then the 
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fragment's rotational velocity will scale as v TOt oc R^ 1 (where R is 
the size of the fragment), and since its free-fall velocity will scale as 
vg oc i? -1 / 2 , this implies that A oc i? -1 / 2 ; in other words, the fragment 
will spin faster as it contracts. The fragment will become centrifugally 
supported once A = 1, and so the contraction of the fragment will stop 
once it reaches a size 

R = \%Ro (31) 

where An, is the initial spin parameter of the fragment, and Ro is its 
initial size. 

Layzer (1963) estimated the initial spin parameter to be of order 
unity, and therefore argued that any fragments that formed would be 
barely distinct from the background gas, and would inevitably collide 
and coalesce before the end of the protogalactic collapse. Hunter (1964), 
however, disagreed and argued that the mutual gravitational torques 
between the fragments would not affect their spins. He noted that in 
a freely-falling, isothermal, inviscid gas, Kelvin's circulation theorem 
would apply, and so therefore the vorticity of a fragment would change 
only due to its expansion or contraction; the gravitational field would 
have no direct effect. In this analysis, Equation (31) still applies, but Ao 
depends purely on the initial vorticity of the gas forming the fragment, 
and is substantially less than one. 

Regardless of the correct value of Ao, it is clear from this analysis 
that at some point the fragments will become centrifugally supported, 
provided that they conserve angular momentum. Once this occurs, the 
assumption of free-fall collapse made by Hoyle, Hunter and many oth- 
ers is clearly no longer appropriate. Moreover, centrifugally supported 
fragments will tend to flatten into rotating disks, which have different 
stability properties from those of collapsing spheres, as discussed at 
some length in Larson (1985). Specifically, one can define a parameter 

<? = ^ < 32 > 

where k is the epicyclic frequency of the disk and ji is its surface den- 
sity, such that disks with Q > Qcrit are stable against gravitational 
fragmentation 3 . The value of Q C rit depends to some extent on the 
equation of state of the gas, but typically Q CI it ~ 0.5-0.6 (Goldreich 
and Lynden-Bell, 1965; Larson, 1985). Furthermore, even if the first 
generation of centrifugally supported fragments remain unstable (i.e. if 
they have Q < Q CI it), subsequent generations will be more highly sta- 
bilized (Larson 1984, 1985), provided that the evolution is isothermal. 



3 Note that this Q is analogous to the Toomre (1964) stability parameter for a 
stellar disk. 
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Therefore, rather than the successive fragmentation envisioned in the 
hierarchical fragmentation scenario, one may instead that find only one 
or two generations of fragments form, with rotation quickly acting to 
suppress fragmentation on smaller scales. 

Another serious critique of the hierarchical fragmentation scenario 
was put forward by Tohline (1980). He pointed out that there is a funda- 
mental difference between the growth of perturbations in a pressure- free 
collapse, such as that considered by Hunter (1962), and in a pressurized 
collapse. In a pressure-free collapse, overdensities are gravitationally 
unstable on all scales, and begin to grow immediately. In a pressurized 
collapse, on the other hand, only overdensities with masses larger than 
the initial Jeans mass can grow to begin with. Although smaller over- 
densities will subsequently become unstable and begin to grow as Mj 
decreases during the protogalactic collapse, the onset of their growth 
is delayed. This effect is particularly pronounced if the initial Jeans 
mass is close to the mass of the protogalaxy, as in this case small 
perturbations will be unable to grow until after the protogalaxy has 
already collapsed by a significant amount. 

Based on this, Tohline argued that one of the fundamental assump- 
tions of the hierarchical fragmentation scenario - the rapid fragmenta- 
tion of the gas down to scales of order of the Jeans mass - is not correct, 
since Mj may vary faster than the gas can respond. In particular, he 
argued that the minimum mass of a fragment at the time that the gas 
becomes optically thick will not be equal to the Jeans mass at that 
moment, since overdensities with M ~ Mj will only just have become 
unstable, and will not yet have had the chance to grow. Instead, the 
minimum fragment mass will correspond to the Jeans mass at some 
earlier time, and will therefore be significantly larger than is predicted 
by the models discussed in the previous section. 

Another way to consider this issue is to examine the dispersion re- 
lation arising from the perturbation analysis. For the classical Jeans 
analysis of plane wave density perturbations in an infinite uniform 
medium, one obtains 

to 2 = c 2 k 2 - 4ttG P , (33) 

and the same relation holds for spherically symmetric density per- 
turbations of the form r _1 sin kr (Larson, 1985). In the pressure-free 
case, c s = 0, and the growth rate of perturbations is scale-free. On the 
other hand, in the pressurized case, c s is non-zero and the growth rate 
increases with decreasing wavenumber, reaching a maximum for k = 0. 
In other words, large-scale perturbations grow faster than small-scale 
perturbations, suggesting that for the latter to win out and to cause 
fragmentation to occur, they must start with substantially higher den- 
sities. Although Equation (33) applies directly only to a rather idealized 
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protogalaxy, it is reasonable to expect to find similar behaviour in the 
more general case. 

A final problem with the hierarchical fragmentation scenario is that 
it assumes that the gas is initially uniform (which implies that t& is the 
same everywhere), and remains so during the collapse. However, this is 
unlikely to be the case. Even if the gas begins its collapse from a uniform 
state, it will tend to become centrally concentrated during the course of 
the collapse (Bodenheimer and Sweigart, 1968; Larson, 1973; Tohline, 
1982) as a pressure gradient builds up between the centre and the edge 
of the protogalaxy. This is significant, because a centrally concentrated 
gas cloud is far more stable than a uniform cloud with regard to the 
growth of small self-gravitating perturbations, as a number of authors 
have demonstrated (Arny, 1966; McNally and Settle, 1980; Silk and 
Suto, 1988; Lacey, 1989). Although many of these analyses assume a 
high degree of symmetry, which makes unclear the extent to which the 
results will apply in more realistic models of collapse, they do provide 
another indication that fragmentation in real protogalaxies will be far 
less effective than the hierarchical fragmentation scenario assumes. 

3.5. Other forms of fragmentation 

In view of the doubts raised in the previous section concerning the 
effectiveness of fragmentation driven purely by gravity, it is worthwhile 
to spend some time examining two other processes which may bring 
about fragmentation of the protogalactic gas: thermal instability and 
supersonic turbulence. 

3.5.1. Thermal instability 

So far we have considered the gas pressure purely as a stabilizing force, 
resisting the action of gravity. However, if the protogalactic gas is ther- 
mally unstable, then the pressure can itself drive fragmentation within 
the gas. 

Thermal instability occurs when small perturbations to the density 
and/or temperature of a region of gas cause the subsequent temper- 
ature evolution of that region to differ significantly from that of the 
unperturbed gas. The case of interest to us is the one in which the per- 
turbations cause accelerated cooling of the perturbed region. This was 
first investigated by Parker (1953), who derived an instability criterion 



where C = (A— r) j p is the specific heat loss function, under the assump- 
tion that the perturbations were isochoric (i.e. that the gas maintained 




(34) 
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a constant density). Field (1965) pointed out that in this case, pressure 
gradients would develop in the gas, and that thermal instability could 
therefore drive dynamical flows. He also argued that on small scales, the 
gas would rapidly respond to any change in temperature by changing 
its density, so as to keep its pressure constant. In other words, small 
perturbations would evolve isobarically rather than isochorically. Field 
(1965) derived an instability criterion for the isobaric case 



where po and To are the unperturbed density and temperature. This is 
the correct criterion for perturbations with wavelengths Af < A < A c , 
where A c = c s i coo i is the distance travelled by a sound wave in a single 
cooling time in the unperturbed medium, and where Af is the Field 
length, given by 



where k is the coefficient of thermal conduction. For perturbations with 
A > A c , Parker's isochoric criterion applies, while for perturbations with 
A < Af, thermal conduction completely suppresses the instability. 

If we apply these criteria to primordial gas, while keeping the chem- 
ical composition of the gas fixed, then we find that the gas is always 
thermally stable. However, if we allow the chemical composition of the 
gas to vary as we perturb the temperature and density, then various 
chemo "-thermal instabilities become possible. The simplest of these oc- 
curs in very hot gas with 10 5 < T < 10 7 K. Within this temperature 
regime, cooling via He II line emission increases sharply with decreasing 
temperature as helium recombines from He in to He n , leading to a 
pronounced thermal instability. This instability has been investigated 
in the context of galaxy formation by Murray and Lin (1990, 1996; see 
also Lin and Murray, 1992, 2000) but seems unlikely to play a significant 
role in the evolution of the first protogalaxies, as the protogalactic gas 
never becomes hot enough to trigger it. 

At lower temperatures, various instabilities associated with the for- 
mation and dissociation of H2 may occur. The first of these was identi- 
fied by Sabano and Yoshii (1977) and analyzed in more detail by Yoshii 
and Sabano (1979). It occurs in gas with a temperature in the range 
2000 < T < 4000 K (with a slight dependence on density) and is caused 
by the collisional dissociation of H2. Within this temperature range, 
the equilibrium abundance of H2 is a sensitive function of temperature, 
due to the strong temperature dependence of the collisional dissociation 
rate (reaction R9). Therefore, gas which has a slightly lower tempera- 
ture than its surroundings will have a higher H2 abundance, and hence 




(35) 




(36) 
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a higher cooling rate. If this leads to a further drop in temperature and 
increase in H2 abundance, then an instability results. Gas cooler than 
about 2000 K does not suffer from the instability because the collisional 
dissociation rate becomes too small to significantly affect the H2 abun- 
dance, which therefore loses its strong temperature dependence, while 
in gas hotter than about 4000 K, the H2 abundance becomes too small 
to provide effective cooling. 

Silk (1983) identified a related instability that appears once the 
three-body formation of H2 becomes effective and occurs for a similar 
reason: a small decrease in the temperature and the associated increase 
in the density lead to an enhanced H2 abundance and higher cooling 
rate, which further perturb the temperature. This instability vanishes 
if the gas becomes fully molecular or becomes optically thick to H2 line 
emission. 

Finally, another potential instability has recently been identified by 
Ripamonti and Abel (2004). This one occurs at very high densities 
(10 14 < n < 10 15 cm -3 ) and is due to a combination of the onset of 
collisionally-induced emission from the H2 (which is highly sensitive to 
the gas density and which quickly comes to dominate the H2 cooling 
rate), and a renewed phase of collisional dissociation within the dense 
gas (which at slightly lower densities is fully molecular). However, the 
resulting instability is very sensitive to the temperature of the g 
can be seen clearly from figures 7 and 8 of Ripamonti and Abel (2004), 
and it remains to be seen whether this instability will actually occur in 
practice. 

Since thermal instability is capable of creating dense structure in 
the gas on all scales larger than the Field length, we might expect it 
to profoundly influence the ability of the gas to fragment. However, 
in practice, it is unlikely to be of great importance in primordial gas. 
There are two main reasons for this. First, the H2-related instabilities 
discussed above all operate within a fairly restricted range of tempera- 
tures. This significantly limits the size of the temperature contrasts that 
can be created, which in turn limits the size of the resulting density con- 
trasts, which can therefore be disrupted more easily by other processes 
such as turbulence (Abel, Bryan, and Norman, 2002). Second, thermal 
instabilities grow on the cooling timescale, t coo \, and will therefore cause 
significantly restructuring of the gas only when i coo i <C idyn- However, 
Omukai and Yoshii (2003) and Ripamonti and Abel (2004) demonstrate 
that thermally unstable protogalactic gas typically has t coo \ ~ £dyn> 
implying that the instabilities do not have sufficient time to grow. This 
conclusion is supported by the results of Abel, Bryan, and Norman 
(2002), whose simulation includes gas in the thermally unstable regime 
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associated with three-body H2 formation, but who find no indication 
that this instability has any significant dynamical effect. 

3.5.2. Supersonic turbulence 

Another way to create dense structures in the gas without the assis- 
tance of gravity is by compressing it in shocks. In particular, if the 
velocity field of the gas is turbulent and supersonic, then large density 
enhancements can be created as the gas is repeatedly shocked. This is a 
process that has received considerable attention in a Galactic context, 
since there is substantial observational evidence for the existence of 
supersonic turbulence in interstellar gas on all scales larger than a 
few tenths of a parsec. Rather than attempting to summarize all of 
this material here - a hopeless task - I refer the reader to the recent 
comprehensive reviews of Mac Low and Klessen (2004), Elmegreen 
and Scalo (2004) and Scalo and Elmegreen (2004), and restrict my 
discussion to a few points of particular relevance to primordial star 
formation. 

Simulations of supersonic turbulence in self-gravitating, isothermal 
gas have been performed by a number of groups (Ostriker, Gammie, and 
Stone, 1999; Klessen, Heitsch, and Mac Low, 2000; Heitsch, Mac Low, 
and Klessen, 2001; Gammie et al., 2003; Bate, Bonnell, and Bromm, 
2003; Li et al, 2004). The gas in these simulations rapidly develops a 
highly inhomogeneous structure, with a density probability distribution 
function (PDF) which is approximately log-normal. Dense, gravita- 
tionally bound cores form in the highest density regions, with a mass 
spectrum that also appears to be log-normal and which resembles the 
stellar IMF. The efficiency with which cores form depends upon the 
properties of the turbulence, and is lower in models with more power 
on smaller spatial scales (Klessen, Heitsch, and Mac Low, 2000), but 
it appears to be very difficult to completely suppress fragmentation: 
this would require strong turbulence on very small scales, which would 
rapidly decay away (Stone, Ostriker, and Gammie, 1998; Mac Low, 
1999) unless driven by some form of energy input on those scales. 

A series of attempts have been made to relate the stellar IMF di- 
rectly to the statistical properties of interstellar turbulence (Larson, 
1981; Fleck, 1982; Elmegreen, 1993; Padoan, 1995; Padoan, Nordlund, 
and Jones, 1997; Myers, 2000; Padoan and Nordlund, 2002), which 
would allow the IMF to be predicted in environments such as early 
protogalaxies for which no direct observational determinations exist. 
However, none of these attempts have met with widespread acceptance, 
and research in this area is still ongoing. 

Extension of these results to primordial gas is further complicated 
by the fact that most studies of turbulent fragmentation assume an 
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isothermal equation of state. This is a reasonable approximation for 
gas in local molecular clouds, since its cooling time is very short, but it 
is unlikely to be appropriate for primordial gas with i coo i ~ idyn- Since 
there are indications that relatively small changes in the equation of 
state can have a large effect on both the shape of the density PDF 
(Passot and Vazquez-Semadeni, 1998) and on the numbers and masses 
of self-gravitating cores that form (Li, Klessen, and Mac Low, 2003), 
a straightforward extrapolation from the isothermal results appears 
unwise. Work in this area is ongoing. 

3.6. Numerical simulations 

As the previous sections hopefully make clear, a number of different fac- 
tors influence the ability of the protogalactic gas to fragment. While we 
can gain considerable insight into the physics of the individual processes 
through the use of simple analytical models, for a proper understanding 
of how the various different processes interact within a real protogalaxy 
we are currently forced to turn to numerical simulations. 

3.6.1. Simulations of local star formation 

Before discussing the results of simulations designed specifically to 
study protogalactic fragmentation and primordial star formation, it 
seems worthwhile to examine what we can learn from simulations de- 
signed primarily to study local star formation. 

One important thing that we have learnt from this kind of simulation 
is the vital importance, when studying gravitational collapse and frag- 
mentation, of resolving the Jeans length throughout the simulation. 
This was convincingly demonstrated by Truelove et al. (1997), who 
showed that if this criterion is not met, then completely artificial frag- 
mentation of the gas can result. This implies that simulations that fail 
to meet this criterion cannot be used to make meaningful predictions 
about gravitationally-driven fragmentation. Although Truelove et al. 
(1997) restricted their attention to grid-based simulations, it has been 
shown that SPH simulations suffer from a very similar problem (Bate 
and Burkert, 1997; Whitworth, 1998). 

Also of interest are the results of a set of simulations performed 
by Tsuribe and Inutsuka (1999ab). They used high resolution SPH 
simulations to study the isothermal collapse of a set of uniform spherical 
clouds with varying ratios of thermal to gravitational energy, parame- 
terized by 



a = 



2GM ' 



(37) 
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where Rq is the initial radius of the cloud, and of rotational to gravi- 
tational energy, parameterized by 

where is the initial angular velocity Previous analytical and numeri- 
cal work (Tohline, 1981; Miyama, Hayashi, and Narita, 1984) suggested 
that such a cloud would collapse to a disk, with a flatness that de- 
pended on the product ao/3o> and that for «oA) < 0.12, this disk would 
subsequently fragment. However, this criterion is clearly not correct 
when (3q is very small, as it predicts that a cloud with ao > 1 should 
collapse and fragment, whereas in reality such a cloud would have a 
mass smaller than the Jeans mass, and would not collapse. Moreover, 
the analytical derivation of this criterion assumes that the collapse is 
homologous, while in reality collapsing clouds would tend to become 
centrally concentrated. 

Tsuribe and Inutsuka find three possible outcomes for their simu- 
lated clouds. When ao and (5q are both small, the cloud collapses to 
a thin disk and fragments, much as was previously envisaged. As ao 
increases, however, the thickness of the disk also increases, and once its 
flatness - defined simply as the ratio of its radius to its scale height - 
falls below a value of approximately 47r, the disk no longer fragments. 
Further increases in ao lead to a final state that increasingly comes to 
resemble the so-called Larson-Penston similarity solution. 4 Finally, for 
sufficiently large ao, collapse is entirely suppressed. 

In the low (3q limit, the boundary between fragmentation and self- 
similar collapse occurs for ao ~ 0.5; with increasing (3q, the boundary 
moves to smaller ao, and is given approximately by ao = 0.55 — (3q for 
< fa < 0.3. 

Although the initial conditions for these simulations are highly ide- 
alized, they do provide strong support for some of the criticisms of 
the hierarchical fragmentation scenario discussed in section 3.1, and 
demonstrate that gravitational fragmentation does not appear to be 



4 The Larson-Penston solution is an asymptotic similarity solution for the isother- 
mal collapse of a sphere, independently derived by Larson (1969) and Penston (1969), 
which describes the collapse at late times and/or at small distances from the centre, 
when the influence of the boundary conditions has become negligible. It can be 
derived numerically from the governing equations of the flow if one assumes that 
the flow is smooth (i.e. that there are no shocks) and that the central velocity is 
zero; this derivation can also be generalized to the case of a polytropic equation of 
state P — Kp 1 (Larson, 1969). Although other similarity solutions are possible (Shu, 
1977; Hunter, 1977; Whitworth and Summers, 1985), the Larson-Penston solution 
appears to provide the best fit to the results of numerical simulations of isothermal 
spherical collapse (Foster and Chevalier, 1993). 
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as effective as originally anticipated, particularly in clouds with large 
values of a®. 

Tsuribe and Inutsuka (2001) extended this analysis, in a more lim- 
ited fashion, to the case of protogalactic collapse. They performed three 
large SPH simulations of the collapse of uniform, spherical protogalax- 
ies, all with the same initial temperature (To = 150 K) and rotation 
parameter (/3o = 0.25), but with differing masses: M = 10 6 , 3 x 10 6 , 
10 7 M , implying that cto = 0.18,0.09 and 0.04 respectively. Rather 
than assuming that the gas remains isothermal, Tsuribe and Inutsuka 
follow its chemical and thermal evolution during the collapse. Based 
on the previous isothermal results, one would expect fragmentation to 
occur in all three simulations, but in fact fragmentation does not occur 
in the 10 6 M protogalaxy, which instead simply forms a single dense 
central core. This difference from the isothermal case is presumably due 
to the fact that cooling is less efficient in these more realistic model pro- 
togalaxies, which therefore evolve as if they had larger values of a®. If 
this interpretation is correct, then it suggests that we can treat Tsuribe 
and Inutsuka's isothermal fragmentation criterion as a necessary, but 
not sufficient, criterion for fragmentation within primordial gas. As we 
shall see below, this interpretation is consistent with the results of more 
detailed protogalactic simulations. 

3.6.2. Filamentary collapse 

The propensity of gravitationally collapsing spheres of gas to settle 
into disks even when /3o is small, noted by Tsuribe and Inutsuka and 
by many other authors, is not unexpected, since we have known for a 
long time that small departures from spherical symmetry are steadily 
magnified during free-fall collapse (Lin, Mestel, and Shu, 1965). More- 
over, flattened, disk-like clouds will generally fragment into filamentary 
structures (Miyama, Narita, and Hayashi 1987ab) which only subse- 
quently fragment into clumps. In light of this, a number of authors 
have considered the problem of filamentary collapse in primordial gas. 

A pressure-supported filament (i.e. one which is not collapsing or 
expanding radially) is gravitationally unstable to perturbations along 
the axis of the filament. Moreover, it is possible to show that for an 
isothermal filament, the fastest growing perturbation is the one with 
wavelength A c ~ ttR, where R is the scale radius of the filament, defined 
as: 



(Stodolkiewicz, 1963). A similar result can be derived for a polytropic 
filament (Larson, 1985). However, filaments formed during protogalac- 




(39) 



and where M is the mass per unit length and po is the central density 
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tic collapse are unlikely to start in hydrostatic equilibrium in the radial 
direction. Uehara et al. (1996) argue that in that case, it is necessary 
to follow the radial evolution of the filament until such time as the 
contraction timescale, given by t con = po/po, exceeds the fragmentation 
timescale, if rag ~ (Gpo)~ 1 / 2 ; at this point, the equilibrium analysis can 
be applied to give an estimate of the resulting fragment mass. 

Uehara et al. (1996) used this approach to study the fragmentation 
of primordial filaments with a range of values of M. They adopted an 
initial temperature, density and molecular fraction appropriate to gas 
that had already undergone significant cooling and collapse within a 
protogalaxy, consistent with their assumption that the filaments had 
formed in a fragmenting disk. They followed the subsequent chemi- 
cal and thermal evolution of the filament in some detail, but treated 
the density evolution in an approximate fashion: the scale radius was 
assumed to evolve as 

2C 

£ = -— [M-M C (T)], (40) 

where M c is the mass per unit length of an equilibrium isothermal 
cylinder (Ostriker, 1964) 

A « T > " ^ < 41 > 

in which case the density then follows from Equation (39) . They studied 
a number of different collapses, with values of M ranging from 1-2 M c . 
They found that as M increased, there was a corresponding increase 
in the density at which the filament fragmented, and a consequent 
decrease in the fragment mass Mp, which fell from 100 Mq for M = M c 
down to 2 M for M = 2M C . 

Nakamura and Umemura (1999) reconsidered this problem and im- 
proved on the Uehara et al. (1996) analysis in several important re- 
spects. Most significantly, they replaced the approximate treatment 
of the density evolution used by Uehara et al. with a more accurate 
treatment based on a one-dimensional hydrodynamic simulation of the 
collapse of the filament. They also made use of a more extensive chem- 
ical model and considered a wider range of initial conditions. They 
found that the collapse led to one of two possible outcomes, depending 
on the initial temperature and the value of / = M c /M. In the majority 
of the models, H2 cooling was effective at the start of the simulation 
and remained so until the gas became optically thick at high density 
(n ~ 10 11 cm~ 3 ). This allowed the filaments to collapse dynamically 
to high density, fragmenting only once cooling became ineffective, and 
producing fragments with masses Mp ~ 1-10 M . However, in models 



first.tex; 2/02/2008; 13:33; p. 37 



38 



Simon Glover 



where the initial temperature was low (To = 100 K), H2 cooling was 
not immediately effective, and the initial evolution of the filament was 
adiabatic. In models where the mass per unit length was large (/ > 2), 
collapse could persist until the temperature increased to a point at 
which H2 cooling became effective, following which the evolution of the 
filament would continue much as if it had started with a higher initial 
temperature. If the initial mass per unit length were small, however, 
then collapse would very quickly come to an end, resulting in an equilib- 
rium filament with a relatively low central density (n ~ 10 5 cm~ 3 ) that 
would produce much larger fragments with masses of a few hundred 
M . 

Nakamura and Umemura (2001) further improved the treatment 
of this problem by performing two-dimensional axisymmetric hydro- 
dynamical simulations of filamentary collapse, which allowed them to 
follow the fragmentation numerically, rather than estimating its effects 
analytically. They again examined a wide range of initial conditions, 
although in light of their previous results, they restricted their attention 
to filaments with initial temperatures To > 300 K. Once more, they 
found two possible outcomes. In filaments with a low initial density 
(n < 10 5 cm" 3 ), fragmentation occurred prior to the onset of three- 
body H2 formation and the resulting fragments were large, with masses 
Mp ~ 100 Mq. On the other hand, in filaments with a high initial den- 
sity, fragmentation is delayed until after the gas has become optically 
thick, resulting in much smaller fragments of mass M-p ~ 1-2 M Q . 

The potential role played by HD molecules in filamentary collapse 
has been investigated by Uehara and Inutsuka (2000) and Nakamura 
and Umemura (2002). Uehara and Inutsuka assumed that the filaments 
formed in a shock-bounded sheet, and therefore adopted initial con- 
ditions appropriate to primordial gas which has cooled rapidly from 
temperatures T » 10 4 K. As previously demonstrated by Mac Low 
and Shull (1986) and Shapiro and Kang (1987), hydrogen recombina- 
tion lags behind cooling in such gas, resulting in an elevated fractional 
ionization that allows a substantial H2 fraction (/h 2 ~ 10~ 2 ) to build 
up. Uehara and Inutsuka demonstrate that in these conditions, the 
fractional abundance of HD would be of order 10~ 5 , and showed that 
this amount of HD is enough to cool the filament to 50 K and to keep 
it evolving isothermally at this temperature until the HD lines become 
optically thick. They argued that this allows very low mass fragments 
to form, with masses Mp ~ 0.01-0.1 M . 

Nakamura and Umemura (2002) examined collapse from a wider 
range of initial conditions than Uehara and Inutsuka (2000), and showed 
that HD cooling would only be significant if the initial H2 and HD 
abundances were both large. However, even in this case, they obtained 
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a much larger fragment mass than Uehara and Inutsuka (2000), finding 
Mp ~ IOMq. They ascribe this difference in part to their more detailed 
treatment of optical depth effects, and in part to the fact that Uehara 
and Inutsuka assumed that the minimum fragment mass would equal 
the Jeans mass, rather than the mass contained within the fastest 
growing perturbation, which in this case is about an order of magnitude 
larger. 

Filamentary collapse has also been studied by Flower (2002) and 
Flower and Pineau des Forets (2003), who examined the effects of 
including a magnetic field. Flower (2002) used a dynamical treatment 
similar to that of Uehara et al. (1996) to show that even a relatively 
weak axial magnetic field would soon provide enough pressure to halt 
the collapse, resulting in the formation of massive fragments, with sizes 
ranging from Mp ~ 60 M for an initial field strength of 10~ 9 G up to 
Mp ~ 6000 Mq for an initial field strength of 10~ 7 G. However, Flower 
and Pineau des Forets (2003) subsequently showed that if the effects of 
ambipolar diffusion were also included, then the field would be far less 
effective at slowing the collapse, since the fractional ionization of the 
gas in the filament is too low to keep the field strongly tied to the gas. 

Ultimately, in spite of the attention paid to filamentary collapse, the 
relevance of these results to fragmentation in real protogalaxies remains 
unclear. The main concerns are that all of these simulations assume 
initial conditions that are far more smooth and symmetrical than will 
actually be the case in a real collapse, and that they neglect a number 
of potentially important effects such as rotation and turbulence. 

3.6.3. Three-dimensional simulations of protogalactic collapse 
Relatively few 3D simulations of protogalaxy formation have been per- 
formed to date, and of these the only ones with sufficient dynamical 
range to study fragmentation within the newly formed protogalaxy 
are the SPH simulations of Bromm, Coppi, and Larson (1999, 2002), 
and the adaptive mesh refinement simulations of Abel, Bryan, and 
Norman (2000, 2002), both of which were discussed previously in sec- 
tion 2.5. 

As previously noted, Bromm, Coppi, and Larson study collapse from 
somewhat idealized initial conditions: a single, isolated overdensity, 
in rigid rotation with spin parameter A, and set to collapse at some 
specified redshift z c . They begin their simulations at z = 100, and 
evolve from then until shortly after the protogalaxy has virialized. By 
focusing on a single protogalaxy in this way they are able to achieve 
a high mass resolution. The precise resolution depends on the mass 
of the protogalaxy simulated and the number of SPH particles used, 
but for their fiducial case of a 2 x 10 6 M Q protogalaxy with a baryon 
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fraction of 5%, simulated with 16384 particles, Bromm, Coppi, and 
Larson achieve a mass resolution of approximately 200 M Q . To avoid 
the numerical difficulties that would otherwise force the simulations to 
be halted once the Jeans mass fell below this value (Bate and Burkert, 
1997; Whitworth, 1998), Bromm, Coppi, and Larson make use of a sink 
particle technique (Bate, Bonnell, and Price, 1995). SPH particles with 
densities greater than a threshold value n tn = 10 s cm -3 and which 
are in regions of converging flow (V • v < 0) are removed from the 
simulation, and replaced with one or more sink particles. One sink 
particle is created for each individual collapsing region, with a mass 
and momentum equal to the sum of the masses and momenta of the 
particles that it has replaced. Once created, sink particles continue to 
interact with the surrounding gas particles, and can accrete them if 
they lie within two smoothing lengths and satisfy the criteria above. 
Further details of the algorithm are given in Bromm, Coppi, and Larson 
(2002). 

Bromm, Coppi, and Larson (1999) present results from a single sim- 
ulation of protogalactic collapse. The protogalaxy they simulate has a 
total mass of 2 x 10 6 M Q , a baryon fraction of 5%, a spin parameter 
A = 0.05, and collapses at a redshift z = 30. Following the initial 
sequence of compression, shock and subsequent cooling and settling 
that has already been described in section 2.5, Bromm, Coppi, and 
Larson find that the cooled gas settles into a flattened central disk, 
with a radius of approximately 15 pc and thickness of 2 pc. This disk 
rapidly breaks up into about a dozen dense clumps, with masses rang- 
ing from 200-10 4 M . The gas in the disk has a mean temperature 
of approximately 200 K (although there is considerable scatter about 
this value), and a density of order 10 4 cm -3 . The gas in the clumps is 
somewhat hotter (T ~ 500 K) and substantially denser, with densities 
ranging all the way up to ra tn . Bromm, Coppi, and Larson (1999) find 
no evidence for further fragmentation of the clumps, but since they 
are soon replaced in the simulation by sink particles, they are unable 
to rule it out. However, further fragmentation of the clumps appears 
unlikely: they have relatively large ratios of thermal to gravitational 
energy («o — 0.5 for a typical clump), and so based on the Tsuribe and 
Inutsuka criterion, one would not expect them to fragment. 

Bromm, Coppi, and Larson (2002) discuss the results of a more 
extensive range of simulations, which sample a wider range of initial 
conditions. They find that the thermodynamic behaviour of the gas is 
very similar in each of their simulations. In every case, the gas cools 
rapidly until it reaches a temperature and density at which cooling 
becomes ineffective. In a gas dominated by H2 cooling, this occurs 
at T ~ 200 K and n ~ 10 4 cm" 3 : below this temperature, the H2 
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cooling rate falls off exponentially, while above this density, collisional 
de-excitation of the excited levels of H2 significantly reduces its effec- 
tiveness as a coolant. Since fragmentation typically occurs after the 
gas has reached this state, the fragment masses tend to lie close to the 
Jeans mass corresponding to this temperature and density. 

In contrast, the morphology of the cool gas is far more sensitive 
to the initial conditions of the simulation. In most cases, it is disk- 
like, but the size and visual appearance of the disks alter significantly 
as the spin parameter and the collapse redshift are varied. A notable 
exception is the single low-mass protogalaxy simulated by Bromm, 
Coppi, and Larson (2002), which had a total mass M = 2 x 10 5 M , 
and the usual baryonic fraction of 5%. This had a larger degree of 
pressure support than the more massive protogalaxies, and settled into 
a spheroidal, quasi-hydrostatic equilibrium configuration at a density 
of order 10 2 cm~ 3 . Bromm, Coppi, and Larson found no evidence for 
fragmentation within this protogalaxy, although a single massive dense 
clump did form in its centre, much as in the high «o simulations of 
Tsuribe and Inutsuka (1999a). 

Finally, Bromm, Coppi, and Larson also examined the effects of HD 
cooling and showed that it made very little difference to the outcome 
of the simulation. This appears to be due to the fact that although 
HD becomes the dominant coolant in low temperature gas, it never 
becomes an effective coolant - the cooling time of the gas remains 
significantly longer than the dynamical time, and so the gas does not 
become significantly cooler than it would if the HD were not included 
in the simulation. 

Abel, Bryan, and Norman (2000, 2002) pursed a rather different 
strategy in their study of protogalactic collapse. Rather than simulating 
collapse from a range of different initial conditions, they instead focused 
on simulating a single example in great detail, starting from realistic 
initial conditions within a large simulation volume, and following the 
collapse to higher densities than those reached in Bromm, Coppi, and 
Larson's SPH simulations. On large scales, their results agree with 
those of other simulations of protogalactic collapse, as I have already 
discussed in section 2.5. On smaller scales, they find an accumulation 
of cold gas within the central ten parsecs of the protogalaxy, much as 
Bromm, Coppi, and Larson do. However, the morphology of this gas 
is not at all disk-like - it is more like a slightly flattened spheroid, 
although less symmetric than this description suggests. Abel, Bryan, 
and Norman find no evidence for any fragmentation of this cool gas, 
beyond the formation of a single dense core of mass M ~ 100 M at 
its centre. 
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At the moment that it forms, this central core has a similar tem- 
perature and density to the surrounding cool gas, namely T ~ 200 K 
and n — 10 4 cm~ 3 , giving it a ratio of thermal to gravitational energy 
ao — 0.5. It is also rotating slowly, with [3q = 0.01. These values 
suggest that even if the core were to collapse isothermally, it would 
be unlikely to fragment further. In fact, the collapse of the core is not 
isothermal: instead, the gas rapidly heats up to a temperature of about 
800K. Indeed, it is this sharp rise in temperature, as much as anything, 
that distinguishes the core from its surroundings, as its density profile 
merges smoothly with the surrounding gas. 

This rise in temperature, which is also seen in the cores that form in 
the Bromm, Coppi, and Larson simulations, albeit to a lesser degree, 
is a natural consequence of the reduced efficiency of H2 cooling at high 
densities. Above a critical density of approximately 10 4 cm -3 , the H2 
cooling rate begins to scale as Ah 2 oc n, while the compressional heating 
rate scales as r comp oc n 3 / 2 , and so the latter eventually begins to 
dominate, causing the core temperature to increase. 

In their original simulation, Abel, Bryan, and Norman followed the 
evolution of the core only up to a density of 10 8 cm -3 ; as their chemical 
model did not include 3-body H2 formation, which becomes effective 
at this density, any results from higher densities would have been 
highly inaccurate. In their subsequent simulation, they included the 
three body reactions, allowing them to follow the collapse of the core 
to much higher densities. They were eventually forced to stop at a 
density of 10 13 cm -3 , because the H2 rotational and vibrational lines 
were becoming optically thick, and their assumption of optically thin 
cooling was therefore no longer valid. 

Abel, Bryan, and Norman found no evidence for fragmentation of 
the central core in either of their simulations. In particular, the thermal 
instability associated with H2 formation and discussed in section 3.5.1 
appears to have only a minor effect on the evolution of the core, and 
does not cause it to fragment. Turbulence is also ineffective at driv- 
ing fragmentation, since the collapse is predominantly subsonic, only 
becoming marginally supersonic within the central 10-20 AU near the 
end of the simulation. 

Finally, Abel, Bryan, and Norman show that the core never becomes 
rotationally supported, and that its final rotational velocity is about 
half of the Keplerian orbital velocity t>Kep = 

(GM/r) 1 / 2 . This is one 
of their most surprising results, as it implies that angular momentum 
is being transferred outwards during the evolution of the core. Indeed, 
Abel, Bryan, and Norman are able to demonstrate directly that this is 
occurring (see figure 4a of Abel, Bryan, and Norman, 2002). This result 
inevitably raises the suspicion that it is due to some purely numerical 
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effect, such as numerical shear viscosity (Norman, Wilson, and Barton, 
1980). However, Abel, Bryan, and Norman find that the details of the 
angular momentum transfer are independent of the hydrodynamic algo- 
rithm used and of the spatial resolution (provided that the simulation 
continues to satisfy the Truelove criterion) . This suggests that the effect 
that they find is real, but of course is not yet conclusive; it would 
be extremely useful to be able to reproduce this result with another 
hydrodynamical code, ideally one based on a fundamentally different 
algorithm, such as SPH. 

Abel, Bryan, and Norman (2002) ascribe the angular momentum 
transfer to the action of shocks during the collapse, but this conclusion 
is open to question since, as noted above, the infall is subsonic in most of 
the core. Another possibility is that angular momentum is transferred 
by tidal interactions with external mass concentrations (Larson, 2002). 
Ultimately, to develop an understanding of the physics underlying this 
effect, we are likely to require additional high resolution adaptive mesh 
simulations (Norman, 2003). 

3.6.4. The optically thick phase 

To follow the evolution of the gas beyond the density reached by Abel, 
Bryan, and Norman (2002), it is necessary to solve a radiative transfer 
problem for the optically thick H2 line emission, in order to be able to 
calculate the correct cooling rate. Unfortunately, an exact solution of 
this problem within a three-dimensional hydrodynamical simulation is 
not currently feasible, since it is essentially a seven-dimensional prob- 
lem (three spatial dimensions, two angles plus frequency and time). 
Approximate methods, such as the otvet formalism of Gnedin and 
Abel (2001) will help in the near future, but so far the only simu- 
lations that have been performed of this last stage of protogalactic 
evolution have been forced to assume spherical symmetry, purely for 
reasons of efficiency. This unfortunately renders them mute on such 
topics of interest as whether the efficient outward transfer of angular 
momentum found by Abel, Bryan, and Norman continues at higher 
densities, or whether the core fragments into a close binary or multiple 
system rather than a single star. 

The first detailed simulations of the evolution of the core in the 
optically thick phase were performed by Omukai and Nishi (1998). They 
used an explicit one-dimensional Lagrangian hydrodynamical code to 
simulate the collapse of a small number of model cores with different 
masses and densities. At gas densities below 10 15 cm -3 , they followed 
the chemical evolution of the gas using a simplified chemical model 
based on that of Palla, Salpeter, and Stahler (1983); at higher densities, 
chemical equilibrium was assumed, and the chemical abundances were 



first.tex; 2/02/2008; 13:33; p. 43 



44 



Simon Glover 



obtained from solution of the coupled Saha equations. Cooling from 
both H2 line emission and collision-induced continuum emission was 
included; the latter dominates at very high densities. Omukai and Nishi 
computed the radiative transfer of this emission using the tangent ray 
method (Hummer and Rybicki, 1971), under the assumption that line 
transfer and continuum transfer can be decoupled. To further simplify 
the calculation, they assumed that the diffusion approximation holds in 
regions that are highly opaque in the continuum, and that line cooling 
from this gas was negligible. 

Omukai and Nishi found that after a short initial transient, the 
evolution of each of their model cores was essentially the same, and 
so they presented detailed results for only a single example: a poly- 
tropic core of mass M = 100 M Q and density at the half-mass radius 
n n = 10 6 cm -3 . As it collapsed, this core quickly developed a self- 
similar density profile with p oc r~ 2 2 , and the collapse as a whole was 
well described by a Larson-Penston type similarity solution for a gas 
with an equation of state p = Kp^, where K = 4.2 x 10 11 (in cgs 
units), and 7 = 1.09. The gas in the centre of the collapsing core soon 
became optically thick in the H2 lines, but this did not immediately lead 
to the evolution becoming adiabatic, as enough cooling was possible 
through the optically thin wings of the lines, as well as in the continuum 
via collision-induced emission, to maintain 7 cr T < 4/3 for an extended 
period. Eventually, however, the core temperature became high enough 
to thoroughly dissociate H2, and the evolution became fully adiabatic. 
This occurred for a central density n c = 10 22 cm -3 , and lead to the 
formation of a small hydrostatic core, with mass M = 5 x 10 -3 M , 
at the centre of the flow. This core rapidly became fully ionized, and 
was bounded by an accretion shock at a radius of 2 AU, and it seems 
natural to identify it as a protostar. Unfortunately, Omukai and Nishi 
were unable to follow its subsequent evolution, as the Courant timestep 
became prohibitively small once the core had formed, forcing them to 
terminate their simulation. 

More recently, Ripamonti et al. (2002) also simulated the optically 
thick collapse phase. Their basic approach was similar to that of Omukai 
and Nishi (1998), but with two major improvements: they included a 
term in the momentum equations corresponding to the force exerted 
on the gas by the scattered H2 emission, and they used a more de- 
tailed model for the chemical evolution of the gas and the behaviour 
of the equation of state at very high densities, based on Saumon, 
Chabrier, and van Horn (1995). In addition, they also examined a 
wider range of initial conditions. Despite this, they found essentially 
the same behaviour as Omukai and Nishi. The evolution of the model 
cores was strongly convergent and soon became well described by a 
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Larson-Penston type similarity solution. This self-similarity lasted until 
a small hydrostatic core of mass 3 x 10~ 3 M Q formed at the centre of 
the flow. 

3.7. Summary 

In the introduction, I posed a number of questions concerning the 
evolution of gas within newly formed protogalaxies, namely: does the 
gas fragment? If so, how large are the fragments? And when and why 
does fragmentation stop? 

Much of the work that has been done on primordial star formation 
assumes some version of the hierarchical fragmentation scenario, in 
which fragmentation is highly efficient and is terminated only by the 
transition of the gas from isothermal to adiabatic evolution. In these 
models, the main uncertainties are the cause of this transition - chem- 
ical changes or fragment opacity? - and the temperature and density 
at which it occurs. 

However, as I outlined in section 3.4, there are a number of rea- 
sons to believe that hierarchical fragmentation does not occur in real 
protogalaxies as various effects combine to inhibit fragmentation. This 
conclusion is supported by the results of the simulations of Tsuribe 
and Inutsuka (1999ab, 2001), Bromm, Coppi, and Larson (1999, 2002) 
and Abel, Bryan, and Norman (2000, 2002): in each of these simula- 
tions there is at most a single episode of fragmentation, and no evi- 
dence for any subfragmentation (i.e. fragmentation of the fragments). 
Moreover, in some of these simulations, such as Abel, Bryan, and Nor- 
man (2000, 2002), the use of the word 'fragmentation' to describe the 
evolution of the gas is misleading: the single 'fragment' that forms 
is actually just the central dense core of a more extended density 
distribution. 

Why is it that fragmentation is so inefficient? Inefficient fragmen- 
tation appears to be a natural outcome of quasi-spherical collapse in 
gas with a high ratio of thermal to gravitational energy. During such 
a collapse, the large thermal pressure will create a strongly peaked 
density distribution, even if the gas is initially quite uniform. It will also 
suppress small-scale collapse until the density of the gas has increased 
by a large factor, since the local Jeans mass scales as Mj oc n -1 / 2 
in isothermal collapse. In combination, these effects imply that the 
gas at the centre of the protogalaxy will quickly come to have both 
the smallest Jeans mass and the shortest free-fall time and therefore, 
unless its collapse is delayed or halted in some way, it will continue to 
collapse all the way to protostellar densities before the bulk of the gas 
has had the opportunity to fragment. If this interpretation is correct, 
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it suggests that further fragmentation may occur within, say, the Abel, 
Bryan, and Norman simulations, if they were continued past the point 
at which the first star forms. However, since this first star will exert a 
strong feedback on its surroundings on a short timescale (Omukai and 
Nishi, 1999; Glover and Brand, 2001; Bromm, Yoshida, and Hernquist, 
2003; Whalen, Abel, and Norman, 2003) it appears unlikely that any 
further fragmentation would in fact occur. 

The efficacy of fragmentation could be enhanced by delaying the col- 
lapse of the densest gas, giving the lower density gas more time in which 
to fragment. As Bromm, Coppi, and Larson demonstrate, rotation can 
do this to some extent, but the outward transfer of angular momentum 
identified by Abel, Bryan, and Norman makes it less effective than sim- 
ple estimates would suggest. Strong perturbations arising from thermal 
instability or supersonic turbulence could also boost fragmentation, but 
in practice neither process is particularly effective in primordial gas. 

The few fragments which do form typically have initial masses of 
a few hundred M . This particular mass scale appears to be a conse- 
quence of the thermodynamics of the gas. At densities less than the 
H2 critical density of 10 4 cm" 3 , H2 cooling is efficient, and the gas can 
cool to a minimum temperature of about 200 K. At higher densities, 
H2 cooling becomes less efficient, and the gas heats up. These values of 
density and temperature therefore mark the point at which isothermal 
evolution comes to an end and ~f c g first exceeds one, and so it is not 
surprising that the minimum fragment mass corresponds approximately 
to the value of the Jeans mass at this density and temperature. 

The major uncertainty that remains concerns the behaviour of the 
collapsing gas in the optically thick regime. It may continue to collapse 
quasi-spherically, in which case we would expect a single, low-mass 
protostellar core to eventually form, as in the simulations of Omukai 
and Nishi (1998) and Ripamonti et al. (2002). On the other hand, it may 
form a gravitationally unstable disk, in which case fragmentation into 
a binary or multiple system would appear to be more likely. Resolution 
of this uncertainty awaits the development of an accurate and efficient 
way of treating the thermal evolution of the optically thick gas. 



4. Protostellar accretion and the final stellar mass 

The results of the simulations of Omukai and Nishi (1998) and Ri- 
pamonti et al. (2002) suggest that the initial mass of a primordial 
protostar may be very small, no more than a few thousandths of a 
solar mass. However, this small initial protostar will be surrounded by 
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a large envelope of infalling gas, some fraction of which will inevitably 
be accreted by the protostar. 

If mass loss from the protostar is negligible (which seems to be a 
good approximation even for very massive metal-free stars - see Marigo, 
Chiosi, and Kudritzki, 2003), then the final stellar mass M* is related 
to the initial protostellar mass M pr by 

M*(t) = M pr + /' M(t)dt. (42) 
Jo 

The final mass is therefore determined by the evolution of the mass 
accretion rate M(t) over the lifetime of the star. This in turn is in- 
fluenced both by the properties of the gas surrounding the star - the 
amount of gas available, its temperature and angular momentum, etc. 
- and by the effects of feedback from the star onto the gas, in the form 
of radiation and outflows. 

4.1. Accretion in the absence of feedback 

Protostellar feedback is complicated to model, and so it is easier to 
begin by considering models of protostellar accretion that do not in- 
clude its effects. Since feedback will act to reduce the accretion rate, 
and hence also the final stellar mass, these models allow us to place an 
upper limit on the ultimate mass scale of the first stars. 

One possible approach to determining the accretion rate is to con- 
struct a simplified model for the collapsing protostellar core from which 
an approximation to the true accretion rate can be derived. For in- 
stance, if we assume that the protostellar core is isothermal and spher- 
ically symmetric, then there exists an entire family of similarity solu- 
tions that could potentially be used to describe the collapse (Hunter, 
1977; Whitworth and Summers, 1985), of which the most familiar are 
the Larson-Penston solution (Larson, 1969; Penston, 1969) and the Shu 
solution (Shu, 1977). 

This approach has recently been applied to primordial star formation 
by Tan and McKee (2004). They model the accretion flow as a spherical, 
isentropic polytrope, and derive an accretion rate that is a function of 
three parameters: the entropy parameter K = p/p* 1 , the polytropic 
index 7 P (which, for an isentropic flow, is equal to the adiabatic index 
7), and </>*, a numeric parameter of order unity, which is related to 
the initial conditions of the flow. Tan and McKee normalize these 
parameters based on the numerical results discussed in the previous 
sections, and set 7 P = 1.1, = 1.43 and K = 1.88 x W 12 K' (in cgs 
units), where 
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and where the effective temperature T e ff = P c g / (nk) includes the small 
contribution to the effective pressure made by subsonic turbulence. 
With these values, the accretion rate becomes 5 

, / f \ -0-30 

M = 7.0 x 10" 2 K /3/ [j^J M © y^ 1 - ( 44 ) 

This is plotted in figure 3 for the case of K' = 1. 

Another obvious approach is to simulate the accretion flow numeri- 
cally, but, as previously discussed, an accurate 3D simulation remains 
impractical due to the expense of the radiative transfer calculations. We 
are therefore forced to choose between simulating the radiative transfer 
and the cooling accurately, at the cost of restricting the hydrodynamics 
to one dimension, or simulating the hydrodynamics correctly, at the 
cost of oversimplifying (or simply neglecting) the radiative transfer 
effects. 

Two important examples of the first approach are the simulations of 
Omukai and Nishi (1998) and Ripamonti et al. (2002), discussed at the 
end of section 3. To recapitulate: Omukai and Nishi use a spherically- 
symmetric Lagrangian code to simulate protostellar core formation 
within initially polytropic clouds, and include a thorough treatment of 
radiative transfer within the H2 lines and in the continuum, using the 
tangent ray method (Hummer and Rybicki, 1971). Omukai and Nishi 
find that prior to core formation the flow is well described by a Larson- 
Penston type similarity solution; specifically, the solution corresponding 
to K = 4.2 x 10 11 (cgs) and 7 = 1.09. They are unable to continue 
their simulations after the formation of the protostellar core, as the 
Courant timestep in the central regions becomes prohibitively small. 
However, if the same Larson- Penston type solution were to apply after 
core formation, then the resulting accretion rate would be 

/ 1 \-0.27 

M = 8.3 x 10" 2 [j^J M © yr _1 (45) 

and the stellar mass would grow as 

/ t \ a73 

M* =0.11^— J M . (46) 

5 Strictly speaking, the accretion rate derived by Tan and McKee (2004) is for 
accretion onto both the protostar and its associated accretion disk. Moreover, they 
also allow for the possibility that protostellar feedback may reduce the amount of 
gas reaching the centre of the system. However, for ease of comparison between their 
results and those of the other authors discussed in this section, I have neglected these 
complications for the time being. 
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The simulations of Ripamonti et al. (2002) are similar in design to 
those of Omukai and Nishi (1998), but incorporate several significant 
improvements. Specifically, Ripamonti et al. include: 

(i) A term in the momentum equation corresponding to the radiative 
force per unit mass 



where n u and F v are the opacity and specific energy flux at fre- 
quency v. 

(ii) An improved treatment of chemistry and thermodynamics at very 
high densities (n > 10 21 cm -3 ), based on Saumon, Chabrier, and 
van Horn (1995), that accounts for non-ideal effects such as pressure 
ionization. 

(iii) A 'frozen core' approximation, which keeps the central mass shells 
fixed in space once their infall velocities fall below a specified value 
(v/vg < 10~ 3 ) and their temperatures exceed 5 x 10 4 K. This 
approximation allows the simulations to avoid the worst of the 
Courant timestep limitations, and hence enables them to be con- 
tinued into the period after core formation. 

Prior to core formation, there is good agreement between the results 
of Omukai and Nishi (1998) and Ripamonti et al. (2002), confirming 
that the flow at this initial stage is well described by a Larson-Penston 
type similarity solution. Differences appear, however, once the proto- 
stellar core has formed. For initial conditions comparable to those stud- 
ied in Omukai and Nishi (1998), Ripamonti et al. obtain an accretion 
rate that is approximately 



which is smaller than the Omukai and Nishi rate and falls off more 
rapidly. 

Ripamonti et al. also find that the accretion rate is sensitive to 
the initial conditions of the simulation: clouds with higher initial tem- 
peratures produce cores with larger accretion rates, even though the 
simulations strongly converge at late times. The difference in rates is 
large enough to produce a difference of a factor of a few in M* by the 
end of the simulations, which follow only the first 10 years after core 
formation. At later times, we would expect the difference to be even 
more pronounced. 




(47) 




(48) 
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Turning to the 'hydro dynamical' approach, the best examples are 
the simulation by Abel, Bryan, and Norman (2002), discussed in de- 
tail in the previous section, and the recent work by Bromm and Loeb 
(2004). In both cases, the full hydrodynamical problem is solved, using 
adaptive mesh refinement in the former case, and SPH with particle 
resampling (Kitsionas and Whitworth, 2002; Bromm and Loeb, 2003) 
in the latter. Additionally, chemistry and cooling are followed accu- 
rately up to the point at which opacity effects begin to dominate. At 
this point, the two treatments diverge. Abel, Bryan, and Norman halt 
their simulation once the optical depth at line centre of the main H2 
cooling lines exceeds 10, at which point the maximum gas density is 
approximately 10 13 cm -3 , and the size of the dense core is a few tens of 
AU. They estimate the subsequent accretion rate based on a calculation 
of the accretion timescale, t acc , at the end of their simulation. They take 
i acc to be 



where v r is the radial velocity of the gas, and they assume that the 
stellar mass at time t is simply the mass of gas with i acc < t. The 
resulting inferred accretion rate is plotted in figure 3. 

Bromm and Loeb (2004) are also unable to follow the flow to very 
high densities, since they too neglect opacity effects when calculating 
their H2 cooling rate. However, unlike Abel, Bryan, and Norman (2002), 
they are not forced to halt their simulation once the gas becomes 
optically thick. Instead, they replace the SPH particles representing 
the dense, optically thick gas with sink particles of the type described 
earlier. Sink particle creation is handled much as in Bromm, Coppi, 
and Larson (2002); the main difference is that the density threshold for 
sink creation is much higher, being set at n t h = 10 12 cm~ 3 . Note that 
while the code is capable of creating multiple sink particles, in practice 
only a single sink is required. By using a sink particle, they sacrifice the 
ability to follow the further evolution of the high density gas, and the 
ultimate formation of the protostar, but in return can continue to study 
the gas flow on larger scales for an extended period. On the assumption 
that all of the gas that is accreted by the sink particle will in reality 
be accreted by the protostar, they derive a protostellar accretion rate 
that is approximated by a broken power-law: 



M{r) M(r) 



(49) 



M(r) A7rp(r)r 2 \v r (r) \ ' 





t > 10 3 yr 



(50) 



first.tex; 2/02/2008; 13:33; p. 50 



Formation of the First Stars 



51 




Time (yr) 

Figure 3. The time-dependent accretion rates predicted by various models of proto- 
stellar accretion. Solid line - Omukai and Nishi (1998); dashed line - Ripamonti et 
al. (2002); dotted line - Abel, Bryan, and Norman (2002); dash-dotted line - Bromm 
and Loeb (2004); dash-dot-dot-dotted line - Tan and McKee (2004). In plotting the 
Tan and McKee rate, I have assumed K' = 1. For t < 40 yr, the predicted accretion 
rate of Abel, Bryan, and Norman (2002) depends on the behaviour of gas on scales 
close to or below the resolution limit of their simulation, and is therefore highly 
uncertain. 



Although the two hydrodynamical models agree well at late times (as 
can be seen by a comparison of their predicted accretion rates, which are 
plotted in figure 3), at early times there is considerable disagreement, 
with Bromm and Loeb (2004) predicting a much higher initial accretion 
rate than Abel, Bryan, and Norman (2002). Unfortunately, with only a 
single example of each simulation available, it is not possible to say how 
much of this disagreement can be attributed to the difference in simula- 
tion methods, and how much simply reflects natural variation between 
the accretion rates in different protogalaxies. Further simulations along 
these lines would clearly be valuable. 

An important open question, which has yet to be studied numer- 
ically, is what role angular momentum plays in the final stages of 
accretion onto the protostar. In particular, we would like to know 
whether angular momentum continues to be transported efficiently 
outwards, as it is in the simulations of Abel, Bryan, and Norman (2002), 
or whether it remains approximately constant once the accretion flow 
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becomes supersonic, as assumed by Tan and McKee (2004). This is 
important because in the former case the bulk of the accretion will 
occur directly onto the surface of the protostar, while in the latter case 
accretion will occur primarily through a circumstellar accretion disk. 

By it very nature, this problem is not one that can be tackled us- 
ing one-dimensional simulations; a full three-dimensional treatment is 
called for. However, disk formation, if it occurs, will take place after the 
flow has become optically thick, since the initial disk radius is predicted 
to be only a few AU (Tan and McKee, 2004). Absent a sudden increase 
in computing power sufficient to allow us to treat the coupled radiative 
transfer and hydrodynamics accurately using an algorithm such as that 
outlined in Hayes and Norman (2003), the best approach is probably 
to look for some approximate treatment that succeeds in capturing the 
essential behaviour of the flow, even if this turns out to be somewhat 
inaccurate. One possible approximation is outlined in Ripamonti and 
Abel (2004). They derive an H2 cooling rate for optically thick gas by 
considering the simplified problem of radiation escaping from a spheri- 
cally symmetric, collapsing protostellar core. The resulting cooling rate 
is well fit by 



where uq = 8 x 10 9 cm -3 , (3 = 0.45, and where Lh 2 thm(T) is the H2 
cooling rate in optically thin gas. Ripamonti and Abel demonstrate 
that this simple approximation performs well in comparison to the full 
radiative transfer solution used in Ripamonti et al. (2002); however, 
its accuracy in the three dimensional case is currently unknown. More 
work along these lines is clearly called for if we are to make progress 
on solving this challenging problem. 

Nevertheless, despite the gaps that remain in our understanding of 
primordial accretion flows, one point stands out clearly: the predicted 
accretion rates are very much larger than those inferred for local proto- 
stars, which are typically of the order of 10 -4 -10 -5 M yr _1 for class 
objects (see, e.g. Maret et al, 2002; Beuther et al, 2002), and which de- 
crease significantly as the protostar evolves (Andre, Ward-Thompson, 
and Barsony, 2000). 

The reason for this difference is straightforward. On purely dimen- 
sional grounds, we would expect the time taken to accrete a mass M 
of gas to be of the order of the free-fall timescale for the gas, unless 
some other effect, such as magnetic support or protostellar feedback, 
were to retard the collapse. Therefore, we can write the mean accretion 
rate for the gas as 




(51) 




(52) 
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where we expect the constant of proportionality to be of order unity. 
For gas with a mean density p, we have oc p~ l l 2 , and hence 

M oc Mp l/2 . (53) 

Now, if this mass of gas is close to being in hydrostatic equilibrium, 
which the results of Abel, Bryan, and Norman (2002) show is a reason- 
able approximation for the gas surrounding the protostellar core, then 
M must be of the order of the Jeans mass; if M <C Mj, the gas would 
not be collapsing, while if M > Mj, it would almost certainly have 
fragmented. Therefore, the accretion rate scales as 

M oc Mj p l/2 , (54) 

or, since Mj oc (T 3 /^) 1 / 2 , 

M oc T 3/2 . (55) 

Since the minimum temperature reached by the primordial gas is more 
than an order of magnitude larger than the temperature characteristic 
of local prestellar cores (which is typically of order 10 K; see Ward- 
Thompson, Andre, and Kirk, 2002), we would expect the accretion 
rate to be correspondingly greater, which is precisely what we find in 
the detailed models discussed above. 

In the absence of significant protostellar feedback, these large pre- 
dicted accretion rates will lead to large final stellar masses. This is 
clearly demonstrated in figure 4, where I plot the final stellar mass 
as a function of time for all of the models discussed above. In every 
case, the final stellar mass grows to more than 100M Q in less than 
10 5 yr. Therefore, unless feedback from the protostar can significantly 
reduce the amount of material that the star accretes over its life- 
time, it will inevitably become very massive, and will either end its 
life as a pair-instability supernova (if its final mass lies in the range 
140 < M < 260 M ), or by collapsing directly to form a black-hole 
(Fryer, Woosley, and Heger, 2001; Heger et al, 2003). The possible 
consequences of this are discussed in some detail by Yoshida, Bromm, 
and Hernquist (2004) and I will not discuss them here. 

4.2. Modeling protostellar feedback 

Any potential form of feedback will be powered, directly or indirectly, 
from one of two sources: the energy released by the infalling matter, or 
the energy produced by nuclear burning within the protostar. To model 
the former, we must model gas flow near the surface of the protostar, 
paying particular attention to the properties of the accretion shock and 
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Figure 4- The time-dependent protostellar masses produced by the accretion rates 
shown in figure 3. 

the circumstellar accretion disk. To model the latter, we must model 
the internal structure of the protostar. In practice, since the dominant 
energy source will change over time, from accretion at early times to 
nuclear burning at late times, an ideal model should treat both regions, 
together with as much of the surrounding gas as possible. 

Unfortunately, computational limitations again restrict us to more 
limited models, and we are forced to approximate. The most significant 
approximation that is commonly made is the assumption of spherical 
symmetry. This is a reasonable approximation for the protostar itself, 
provided rotational effects are not significant, but it does not allow 
us to treat any processes involving the accretion disk, and is therefore 
rather limiting. On the other hand, it dramatically reduces the compu- 
tational requirements of the problem, and consequently continues to be 
a widely used approximation. Indeed, the only treatment of primordial 
protostellar structure and feedback of which I am aware that does not 
assume spherical symmetry is the recent work of Tan and McKee (2004) 
and Tan and Blackman (2004), which is discussed in some detail later 
in this section. 

It is also common to further simplify the problem by splitting it 
into two pieces, and considering the evolution of the structure of the 
protostar (which will strongly influence the strength of any feedback) 
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separately from the effect of feedback on the flow. In other words, 
studies of protostellar structure generally assume a constant accretion 
rate, while calculations focused on the effects of feedback on the accre- 
tion flow frequently assume a constant energy source. This separation is 
purely pragmatic; it is easier to study the different processes separately 
before combining them in more realistic coupled models. 

4.2.1. The evolution of protostellar structure 

The evolution in the structure of a primordial protostar as it accretes 
matter from its surroundings was first studied in detail by Stahler, 
Palla, and Salpeter (1986a). Their strategy followed that of Stahler, 
Shu, and Taam (1980ab, 1981), who had previously studied a similar 
problem for the case of a low-mass, population I star. 

They assume that the accretion process can be treated as a series 
of quasi-steady-state accretion flows onto a hydrostatic core, which is 
bounded by a strongly radiating accretion shock. Within the core, the 
standard stellar structure equations are solved, with the assumption 
that deuterium burning is the only source of nuclear energy. Outside 
the core, the treatment depends on the optical depth of the gas. If the 
gas is optically thin to the radiation from the accretion shock, then the 
accretion flow is assumed to be in free-fall. Otherwise, a more detailed 
calculation is made that incorporates the effects of the radiation force 
on the infalling gas. The accretion shock itself is treated as a simple 
discontinuity; no attempt is made to model its structure in any detail. 
Since the thickness of the shock is small compared to the size of the 
core, this should be a good approximation. 

Stahler, Palla, and Salpeter (1986a) begin with an initial core mass 
of 0.01 M , and give the core an arbitrary initial distribution of specific 
entropy: 



where (3 = 7.39 and so is calculated from the adopted central tempera- 
ture (T c = 10 5 K) and density (p c = 0.28 g cm -3 ) using an equation of 
state taken from Eggleton, Faulkner, and Flannery (1973). The outer 
boundary condition is fixed by the accretion rate, which Stahler, Palla, 
and Salpeter take to be constant, with a value 4.41 x 10 -3 M yr _1 . 

Starting from these initial conditions, Stahler, Palla, and Salpeter 
calculated the subsequent evolution of the protostar until the core mass 
reached a value of 10.5 M . They found that the evolution could be 
divided into three qualitatively distinct phases. 

In the first phase, which lasts until the core mass reaches 0.1 M , the 
protostar relaxes from its initial entropy profile into one appropriate for 
the particular choice of accretion rate. This 'decay of transients' phase 




(56) 
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indicates that although the initial conditions are probably incorrect 
in detail, the flow soon loses all memory of them, and therefore any 
inaccuracy at this stage is unlikely to affect the later results. 

Once the initial transients have died away, the protostar enters the 
second phase of its evolution. During this phase, its central temperature 
remains low (T c ~ 10 5 K), resulting in a high interior opacity and 
hence a low interior luminosity. Consequently, the evolution of the core 
during this phase is almost adiabatic; although the core continues to 
gradually contract, this contraction does not lead to any increase in 
the central entropy. Since the postshock entropy increases over time 
due to the increasing strength of the accretion shock (which is itself a 
natural result of the increasing protostellar mass), the core develops an 
off-centre distribution of entropy and temperature. 

The gas surrounding the accretion shock remains optically thick 
throughout this period. This is a direct result of the high accretion 
rate, which produces a highly luminous accretion shock. This produces 
sufficient radiation to partially ionize the preshock gas in the vicinity of 
the shock, creating a structure known as a radiative precursor. The H~ 
opacity of the dense, partially ionized gas in this radiative precursor 
is more than sufficient to make it optically thick. Stahler, Palla, and 
Salpeter show that the core radius during this period evolves as 

/ M \ 0,27 

r * = 4SA {m^J i^ 1 R • (r,7) 

where Mq = 4.41 x 10~ 3 M Q yr -1 , while the photospheric radius evolves 
as 

*-"(£n£F»* 

so R p > R* throughout. The strong H~ opacity also keeps the photo- 
spheric temperature low (T p ~ 5000 K), which prevents the protostar 
from being able to ionize material outside of its photosphere. 

This near-adiabatic accretion phase comes to an end once the cool- 
ing time of the core, given approximately by the Kelvin-Helmholtz 
timescale 

GM * 

*KH = -^-j-, (59) 

becomes comparable to the accretion timescale t acc = M*/M. This 
occurs for a core mass M ~ 1M Q , and results in the core entering a 
phase of homologous collapse, while energy and entropy are transferred 
outwards in the form of a 'luminosity wave'. The radial position of the 
luminosity peak moves outwards towards the accretion shock, reaching 
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it at about the time that the core mass has reached 8 M . This results 
in a rapid swelling of the outermost layers, which weakens the accre- 
tion shock and leads to it becoming optically thin. Stahler, Palla, and 
Salpeter terminate their simulation shortly afterwards, once the core 
mass has reached 10.5 M . 

Although Stahler, Palla, and Salpeter include deuterium burning as 
a possible energy source, in practice they find that it plays no role at this 
stage of the protostar's evolution, as its central temperature remains 
too low to ignite deuterium. However, since the central temperature 
and density are both rising sharply at the end of the simulation as 
the central regions of the core collapse homologously, it is reasonable 
to assume that deuterium ignition will soon take place. Stahler, Palla, 
and Salpeter (1986b) study the onset of deuterium burning and the 
later onset of hydrogen burning in a subsequent simulation of the pre- 
main sequence evolution of a 5 M primordial protostar. Their initial 
conditions are taken from Stahler, Palla, and Salpeter (1986a), but the 
accretion rate is now set to zero. Stahler, Palla, and Salpeter (1986b) 
find that deuterium ignites approximately 6000yr after the beginning of 
their simulation, with hydrogen ignition following after 2 x 10 5 yr. The 
protostar eventually reaches the zero-age main sequence approximately 
10 6 yr into the simulation. 

An alternative treatment of these later stages of evolution that does 
not assume a negligible accretion rate is that of Omukai and Palla 
(2001). They construct their simulation in the same way as Stahler, 
Palla, and Salpeter (1986a) and assume the same constant rate. The 
only significant technical differences between the two simulations are 
that Omukai and Palla use zero metallicity opacities from Lenzuni, 
Chernoff, and Salpeter (1991) and Iglesias and Rogers (1996) in place 
of the older values used by Stahler, Palla, and Salpeter, and that they 
begin their simulation at the start of the optically thin phase, when 
the core mass has already reached M = 8M . However, unlike Stahler, 
Palla, and Salpeter, they do not halt their simulation once the core 
mass reaches 10.5 M Q ; instead, they continue until well after hydrogen 
ignition. 

They find that the period of optically thin evolution identified by 
Stahler, Palla, and Salpeter (1986a) lasts for only a short time; the 
core reaches a maximum radius of 220 R for a core mass of 11.5 M , 
but shortly afterwards begins a sustained process of contraction. The 
radiative precursor reappears once the core mass reaches 12.4 M and 
persists for the remainder of the simulation. As in the previous optically 
thick phase, H~ opacity keeps the photospheric temperature low. 

Within the core, deuterium burning begins once the mass of the core 
reaches 12 M (corresponding to a time t = 1000 yr after the beginning 
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of the simulation, given the assumed accretion rate), and is complete 
by the time the mass has reached 30 M© (corresponding to t = 5000yr). 
It does not contribute significantly to the protostellar luminosity, and 
has little effect on the structure of the protostar. 

Hydrogen ignition follows once the core mass reaches 80 M Q (corre- 
sponding to t = 1.6 x 10 4 yr). At the same time, the internal luminosity 
nears the Eddington value 



triggering a second phase of expansion. The outer layer of the core 
moves out from 10 R© to 100 R©, although it remains well within the 
photosphere, which has a radius of approximately 1000 R Q at this time. 
As the core expands, the accretion luminosity falls and the radiation 
force on the outer layers of the core declines. It soon becomes too small 
to maintain the expansion, and so the core begins to contract rapidly 
for a second time. From this point on, however, nuclear burning makes 
a substantial and increasing contribution to the total protostellar lumi- 
nosity, which soon reaches Z/Edd for a second time. This triggers another 
phase of radiation-driven expansion, which this time is strong enough 
to halt the accretion. This occurs once the core mass has reached 
M ~ 300 M Q , and Omukai and Palla terminate their simulation at 
this point. 

In order to assess the dependence of this result on the adopted 
accretion rate, Omukai and Palla (2003) performed a similar analysis 
for a range of different values of M. They defined a fiducial accretion 
rate M n d = 4.41 x 10~ 3 M yr" 1 corresponding to the value adopted 
by Stahler, Palla, and Salpeter (1986a) and Omukai and Palla (2001), 
and studied models with rates M = (0.25,0.5, 1.0,2.0) x M n d, a s well 
as a model using the time-dependent accretion rate predicted by Abel, 
Bryan, and Norman (2002). 

The earliest stages of protostellar evolution are qualitatively the 
same in all of these models: we see the same sequence of adiabatic 
growth, propagation of a luminosity wave that triggers expansion of 
the outer layers, followed by rapid contraction. There are quantitative 
differences; for instance, protostars with a larger M have a larger radius 
at a given mass. However, significant differences in behaviour do not 
become apparent until the end of the rapid contraction phase. In the 
fiducial case, Omukai and Palla (2003) confirm their previous result: 
they find two episodes of radiation-driven expansion, the second of 
which is strong enough to terminate accretion onto the protostar. In the 



-^Edd 
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M = 2Mfid case, however, they find that the initial phase of expansion 
is strong enough to halt the accretion, thanks to the larger accretion 
luminosity associated with the larger accretion rate. Consequently, the 
final protostellar mass is smaller, being approximately 90 M . 

In the two models with M < M nc i, however, the outcome is rather 
different. Core contraction comes to an end shortly after the onset of 
hydrogen burning, but there is no subsequent phase of radiation-driven 
expansion, as the protostellar luminosity is never more than 70% of 
Ledd- Instead, the core relaxes quickly onto the zero-age main sequence 
(ZAMS), continuing to accrete all the while. 

Omukai and Palla (2003) show that there is a critical accretion 
rate Merit that separates these two outcomes. Protostars with M < 
M cr i t can reach the zero-age main sequence while still accreting, and 
can therefore grow to extremely large masses, while protostars with 
M > Merit will undergo radiation-driven expansion before reaching 
the ZAMS, and will therefore have smaller masses. To evaluate M cr it, 
Omukai and Palla equate the total luminosity of a protostar that has 
just reached the ZAMS with the Eddington luminosity: 

^Edd = ^zams H 5^ — — , (62) 

-KZAMS 

where the second term on the right-hand side represents the accretion 
luminosity. This equation can be rewritten as 

■ AttcRzams f, L ZAMS \ 

Merit = 1 ; , (63) 

Kcs \ -^Edd / 

where K es is the electron scattering opacity. Evaluating this, we find 
that Merit — 4x lO -3 M0yr~ 1 , coincidentally close to Mm. In principle, 
Mrit has a dependence on the mass of the protostar, but in practice 
this dependence is weak and may be neglected. 

The final model that Omukai and Palla (2003) consider is one with 
a time-dependent accretion rate taken from Abel, Bryan, and Norman 
(2002). In this model, the accretion rate is initially much larger than 
Mrit, but decreases with time, and falls below Mrit when the core mass 
reaches 95 M Q . This initial behaviour of this model is very similar to 
that of the model with M = Mgd, but the two diverge during the rapid 
contraction phase; the time dependent model undergoes a very brief 
period of radiation-driven expansion, but thereafter re-contracts, and 
relaxes onto the zero-age main sequence, following which its evolution 
is indistinguishable from that of the other models with M < M CT it- 

An alternative view of the evolution of protostellar structure is pre- 
sented by Tan and McKee (2004). In contrast to previous authors, they 
do not assume spherical symmetry, allowing them to treat the case of 
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accretion via a circumstellar disk. Tan and McKee fix the size of the 
disk by assuming that angular momentum is conserved by gas within 
the sonic point of the accretion flow. This allows them to write the disk 
radius as 

(64) 

where r sp is the radius of the sonic point, M sp is the mass interior to the 
sonic point, m*<j is the mass interior to (i.e. the mass of the protostar 
plus the disk) and ,/kep is the ratio of the rotational velocity of the gas 
to the Keplerian orbital velocity ^Kep = (GM/r) 1//2 , evaluated at the 
sonic point. Tan and McKee fix r sp using their analytic accretion flow 
solution, discussed in the previous section, and adopt fxep = 0.5, based 
on the results of Abel, Bryan, and Norman (2002). They show that in 
this case, the accretion disk radius becomes 

9/7 

K'- 10 ' 7 AU, (65) 

where K' is given by Equation (43), and where e*d is the fraction of 
the infalling gas that reaches the disk or the star. In the absence of 
protostellar feedback, e*d = 1- 

Provided that 3> r*, which will generally be the case, the bulk of 
the gas will accrete first onto the disk and only later onto the protostar. 
Therefore, the accretion rate onto the protostar, and hence its struc- 
ture, will be determined in large part by the behaviour of the disk. To 
determine the disk structure, Tan and McKee (2004) use the standard 
theory of steady, thin viscous accretion disks (as outlined in Shakura 
and Sunyaev, 1973 or Frank, King, and Raine, 1995), with a spatially 
constant viscosity parameter a. The dominant source of this viscos- 
ity and the appropriate value for a remain uncertain, much as they 
do in the analogous situation in present-day star formation. Possible 
sources of viscosity include gravitational instabilities within the disk 
(Larson, 1984; Lin and Pringle, 1987; Bodenheimer, 1995; Nomura and 
Mineshige, 2000; Gammie, 2001; Johnson and Gammie, 2003), tidal 
interactions with external mass concentrations (Spruit, 1987; Larson, 
1990; Lin and Papaloizou, 1993; Blondin, 2000; Larson, 2002), and 
turbulence generated by the magnetorotational instability (MRI; see 
Balbus and Hawley, 1991, 1998) 

The last of these will only operate if a sufficiently strong magnetic 
field is present in the disk. Kulsrud et al. (1997) showed that a very 
small seed field could be produced via the Biermann battery mech- 
anism (Biermann, 1950) during the collapse of the protogalactic gas, 
and Tan and Blackman (2004) show that although this field is initially 
too small to drive the MRI, a dynamo process acting in the disk will 
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rapidly amplify the field, which soon becomes strong enough to drive 
the instability. In view of this, Tan and McKee (2004) examine the 
behaviour of a disk with a = 0.01, the appropriate value for a disk 
susceptible to MRI (Balbus and Hawley, 1998). 

Having determined the disk structure and the rate at which gas 
accretes from the disk onto the protostar, Tan and McKee (2004) then 
solve for the evolution of the protostellar structure using a modified 
version of the analytic approach developed by Nakano, Hasegawa, and 
Norman (1995) and Nakano et al. (2000). In this approach, the proto- 
stellar radius is found balancing the rate of accretion of energy with the 
rate of change of the total protostellar energy. The internal structure 
of the protostar is not solved for explicitly; instead, it is approximated 
as a polytrope, with the polytropic index fixed by comparison with the 
results of Stahler, Palla, and Salpeter (1986a) and Omukai and Palla 
(2001). 

Tan and McKee (2004) show that if an accretion disk is not present 
(i.e. if /kop = 0), and if M = M^, then this model successfully repro- 
duces the results of Stahler, Palla, and Salpeter (1986a) and Omukai 
and Palla (2001). They also demonstrate that the presence of a disk 
has a relatively small effect on the evolution of the protostar. The 
protostellar radius tends to be somewhat larger than in the spherical 
accretion case, but the protostar still evolves through the same progres- 
sion of adiabatic growth, terminated by the emergence of a luminosity 
wave, followed by rapid contraction that ends once the protostar reaches 
the zero-age main sequence. The major difference from the spherical 
case is in the behaviour of the photosphere. Because most of the gas 
accretes onto the protostar via the disk, the gas density is significantly 
reduced in regions near the protostar but out of the plane of the disk. 
Consequently, the optical depth of these regions is also significantly 
reduced, with the result that the flow becomes optically thin early in 
its evolution. For example, in the model with faep = 0.5, the photo- 
sphere vanishes once the protostellar mass reaches 1 M and does not 
subsequently reappear. As we will see below, this may have a major 
influence on the effectiveness of radiative feedback from the protostar. 

4.2.2. The effects of feedback 

In order for the protostar to substantially reduce the rate at which 
it accretes, it must be able to transfer a significant amount of energy 
and/or momentum to the infalling gas. A number of possible mecha- 
nisms that accomplish this have been suggested, which fall under two 
broad headings: radiative feedback, where radiation from the protostar 
(or the accretion disk) is responsible for transferring energy and mo- 
mentum directly to the infalling gas, and mechanical feedback, where 
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the protostar transfers energy and momentum to some form of outflow, 
which subsequently transfers it to the infalling material. 

In local star-forming regions, the dominant mechanism is a form 
of radiative feedback: radiation pressure exerted on the infalling dust 
grains by the protostar results in a substantial momentum transfer to 
the gas and prevents massive stars from forming, unless the accretion 
rate is very large (Wolfire and Cassinelli, 1987). In dust-free primordial 
gas, however, this process is clearly inoperative, and we must examine 
other possibilities. 

One obvious possibility is that radiation within the optically thick 
rotational and vibrational lines of H2 may exert sufficient pressure to 
slow or stop the infall. However, this seems unlikely to be the case: 
in their simulation of protostellar formation, Ripamonti et al. (2002) 
compute the total opacity of the H2 lines and show that it is never 
more than 5% of the electron scattering opacity, implying that the 
protostar would have to radiate a total luminosity in the H2 lines that 
was many times larger than the Eddington luminosity for this effect to 
be dynamically significant. 

A more interesting possibility is that the buildup of an H 11 region 
around the protostar may terminate the accretion. This idea was first 
discussed in the context of present day star formation by Larson and 
Starrfield (1971), and was re-examined in the context of primordial 
star formation by Omukai and Inutsuka (2002). The basic mechanism 
is straightforward: as the protostar ionizes the gas, it transfers to it a 
considerable amount of thermal energy. If the H II region can expand to 
a radius at which this thermal energy exceeds the gravitational binding 
energy of the gas, then the ionized gas outside this radius will become 
unbound from the central protostar, and little if any of it will ultimately 
be accreted. To assess the effectiveness of this mechanism, we need to 
answer two basic questions: one, does an H 11 region actually form? And 
two, if an H II region does form, can it expand sufficiently to unbind the 
gas, or will it instead be confined to a small radius by the inflow? 

The answer to the first of these questions will depend on the effective 
temperature of the protostar. An isolated, massive metal-free star on 
the main sequence will have an effective temperature of approximately 
10 5 K (Cojazzi et al, 2000), and will emit a substantial number of 
ionizing photons, so in this case it is clear that an H II region will 
form. On the other hand, in the accretion models of Stahler, Palla, 
and Salpeter (1986a) and Omukai and Palla (2001, 2003) discussed 
in the previous section, the protostar is hidden within a much larger 
photosphere, with an effective temperature of only 6000 K, and so no 
H II region will form until the photosphere vanishes at late times. 
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Regarding the second question, Omukai and Inutsuka (2002) show- 
that if the accretion flow onto the protostar is steady and spherically 
symmetric, then an H II region can expand sufficiently to unbind the 
surrounding gas only if it is powered by a flux of ionizing photons 
that exceeds a critical value, Qcrit- They demonstrate that in order to 
calculate Qcrit correctly, it is necessary to take into account an addi- 
tional type of radiative feedback - the force arising due to radiation 
scattering within the H II region. There are two main components to 
this force. One is due to Thomson scattering, and will be negligible 
until the protostellar luminosity approaches L^dd- The other comes 
from the momentum transfer that occurs during photoionization, and 
for a gas in photoionization equilibrium the force per unit mass is given 
by (Haehnelt, 1995) 

hu ion n c n u+ 

F ia d = «B , (66) 

c p 

where «b is the case B recombination coefficient and hv im ~ 13.6 eV 
is the mean energy of an ionizing photon. This radiative force acts 
to reduce the infall velocity within the Hn region, which leads to an 
increase in the density of the ionized gas, since the assumption of steady 
flow implies that p oc where v is the infall velocity. This increased 
density leads in turn to a higher recombination rate, which limits the 
expansion of the Hn region. Omukai and Inutsuka show that the net 
effect is to make Qcrit very large; they find a value 

*"- 6 - 4xlo52 (lpfe)"(lo^) 2 ' < 67 » 

where i?j n is the inner radius of the H n region. For reasonable values for 
the stellar parameters, this gives a value of Q cr it that is about a hundred 
times larger than the actual ionizing flux, suggesting that even if an 
H II region forms, it will be unable to halt accretion onto the protostar. 
It is worth noting, however, that the evolution of the H n region is likely 
to be very sensitive to the density distribution near the protostar and 
it is not clear that this conclusion will still hold if we relax some of the 
simplifying assumptions made above. It would be instructive to redo 
this calculation using a more realistic dynamical model. 

A final form of radiative feedback that has attracted serious con- 
sideration is the scattering of Lyman-a photons by the infalling gas. 
Within the Hll region, this is of only minor importance, but in the 
surrounding Hi gas, where the optical depth to Lyman-a scattering 
is much higher, it may be far more significant (Braun and Dckcl, 
1989; Bithell, 1990; Haehnelt, 1995). Doroshkevich and Kolesnik (1976) 
argue that the radiation pressure exerted by the Lyman-a photons 
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will rapidly expel all of the Hi gas near the protostar, thereby limiting 
the final stellar mass to 10 M or less. On the other hand, Omukai 
and Inutsuka (2002) contend that their treatment overestimates the 
density of Lyman-a photons and so overestimates the effectiveness of 
this mechanism. More recently, Tan and McKee (2003) have presented 
preliminary results from a calculation of the effects of Lyman-a scatter- 
ing that assumes a rotating, axisymmetric inflow. Their results appear 
to support the Doroshkevich and Kolesnik (1976) picture, although 
they quote a larger mass limit of order 20 M . However, full details of 
their calculations have not yet been published, so it is not possible to 
assess the strength of their argument. 

The effects of feedback in the form of protostellar outflows have been 
less well studied than the various radiative effects discussed above. In 
the case of radiatively driven outflows, such as stellar winds from O 
type stars, this neglect is easy to understand - these outflows are driven 
primarily by the scattering of photons in the resonance lines of metal 
ions, and thus grow substantially weaker as the metallicity declines 
(Kudritzki, 2002). In primordial gas, outflows of this type must rely 
on Thomson scattering and thus will only become significant if the 
protostellar luminosity reaches LeoM) which, as we have seen, will only 
occur if the accretion rate exceeds M C rit- 

Meanwhile, bipolar outflows of the kind that are ubiquitous in lo- 
cal star-forming regions have attracted little study because they are 
widely believed to be hydromagnetic in nature (see, e.g. Matzner and 
McKee, 1999) and protogalactic magnetic fields were thought to be too 
weak to power them. However, in a recent paper, Tan and Blackman 
(2004) have argued that the initial protogalactic magnetic field will 
undergo substantial amplification by a helical dynamo operating in the 
turbulent accretion disk surrounding the protostar, and may therefore 
become strong enough to drive an outflow. For a reasonable choice of 
parameters, they find that an outflow with mechanical luminosity 



can be produced. Although this value is substantially less than the 
radiative luminosity of the protostar, the outflow is far more effective 
than the radiation at transferring momentum to the surrounding gas. 
Tan and Blackman estimate that it will begin to remove a significant 
quantity of gas from the core once the protostellar mass exceeds 20 Mq, 
and that as much as 50% of the core mass may have been removed by 
the time that the protostellar mass reaches 100 M Q . The outflow will 
also alter the density structure of the core and should therefore be taken 
into account when assessing the effectiveness of radiative feedback. 




(68) 
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4.3. Summary 

Although much work remains to be done on developing a detailed un- 
derstanding of primordial protostellar accretion, several key points are 
already clear. First, the mean accretion rate of a primordial protostar 
is much larger than that of its present-day counterparts, simply as a 
result of the higher gas temperature of the protostellar core. Second, a 
large quantity of gas is available to be accreted by the protostar, since 
the very low efficiency of fragmentation means that it does not have to 
compete with a large number of other protostars for the available gas. 
Third, many of the forms of protostellar feedback that serve to limit the 
masses of protostars forming at the present day either do not operate 
in the primordial case, or operate with a reduced effectiveness. At the 
same time, the higher accretion rate implies a larger ram pressure of 
the infalling gas, making the job of halting the accretion harder than 
at the present day. Fourth, those forms of feedback that do seem to be 
effective (Lyman-a scattering, hydromagnetic outflows) do not become 
so until late times, after the protostar has already accreted a substantial 
quantity of gas. 

Taken together, these points strongly suggest that the first stars will 
be very massive. Indeed, if this basic picture is correct, it is difficult to 
see how accretion could be terminated early enough to produce a solar 
mass star, since the predicted accretion rates discussed earlier suggest 
that this mass of gas will build up in only 10-20 yr. 

The major uncertainties that remain are easily summarized: 

(i) Is our basic picture of a single protostar per core correct, or does 
a second stage of fragmentation occur at late times, after the gas 
has become optically thick? 

(ii) Does a dynamically significant disk form? If so, how does it evolve, 
and how does it affect accretion onto the protostar? 

(iii) Is the effective temperature of the protostar 6000 K (as suggested 
by the models of Stahler, Palla, and Salpeter, or Omukai and Palla) 
or 10 5 K (as suggested by the models of Tan and McKee)? 

(iv) Are there other possible forms of feedback that we haven't yet 
considered? 

As with many of the open questions concerning primordial star for- 
mation, resolution of these issues is likely to require detailed numeri- 
cal simulations with an adequate treatment of the effects of radiative 
transfer. 
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5. Conclusion 

The past ten years have seen significant advances in our knowledge 
of many aspects of primordial star formation, from the large scale 
environment in which it occurs to the very small scale structure of the 
first protostellar core. Although a number of issues remain unresolved, 
a consensus now exists on the broad outlines of the process. 

We expect the first stars to form in small, H2-cooled protogalaxies, 
with masses of 10 5 -10 6 M , at redshifts z = 30-40. Fragmentation 
of the gas within these protogalaxies will be inefficient, contrary to 
previous expectations, and in the smallest protogalaxies only a single 
dense core will form, with a mass of a few hundred solar masses. 

This core will collapse without fragmenting further until it becomes 
optically thick. Its subsequent evolution is not entirely certain, but the 
most probable outcome is the formation of a single low-mass proto- 
star near its centre. This protostar will accrete gas rapidly from its 
surroundings, and will soon become very large. Protostellar feedback 
may act to limit the accretion rate at late times, in which case the 
final mass of the star will be similar to that of a Galactic O type star; 
otherwise, the final stellar mass will be limited only by the amount of 
gas available, and will be of the order of a few hundred solar masses. The 
first stars will therefore end their lives either exploding as supernovae, 
or collapsing directly to form black holes. Either way, there should be 
none left alive at the present day. 

The future also holds great promise for the study of primordial star 
formation. New facilities such as alma and jwst will for the first time 
allow us to probe the earliest epochs of star formation, and may allow us 
to test observationally the picture that I have outlined above (although 
the practical difficulties will remain formidable). Meanwhile, further 
increases in computing power will allow us to perform increasingly 
detailed simulations and should soon allow us to fill in many of the 
gaps in our current understanding of the formation of the first stars. 
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